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EXECUTIVE  SUMMARY 


DNAPL  AND  LNAPL  DISTRIBUTIONS  IN  SOILS:  EXPERIMENTAL  AND 

MODELING  RESULTS 

By 

D.S.  Dumford,  D.B.  McWhorter,  F.  Marinelli,  C.D.  Miller, 

A.  Swanson  and  H.  Trantham 


The  Air  Force  Office  of  Scientific  Research  Program  in  Environmental  Quality  supports 
research  that  provides  a  technology  base  for  eliminating  potential  environmental  hazards. 
The  Air  Force  uses  a  variety  of  hazardous  materials  that  can  be  introduced  into  the  soil 
environment  through  spills  or  leaks.  Of  interest  in  this  report  are  the  thousands  of  sites 
contaminated  with  fuels  or  solvents.  These  fluids  are  classified  as  nonaqueous  phase 
liquids  (NAPLs)  because  they  have  low  solubilities  in  water  and  their  fate  and  transport 
in  soils  are  determined  by  the  physics  of  immiscible  fluid  flow  in  porous  media. 
Understanding  immiscible  fluid  flow  in  porous  media,  with  the  goal  of  characterizing  and 
remediating  subsurface  sites  contaminated  with  nonaqueous  phase  liquids,  is  the  focus  of 
this  project. 

This  report  summarizes  three  different  components  of  AFOSR  Project  No.  F49620-94-1- 
0193.  Section  II  of  the  report  provides  examples  of  seemingly  subtle  aspects  of 
immiscible  fluid  behavior  that,  nonetheless,  are  manifested  in  important  large-scale 
effects.  Results  of  one  and  two-dimensional  multiphase  laboratory  experiments  in 
homogeneous  and  layered  soils  are  presented.  Results  show  the  importance  of  hysteresis 
in  the  capillary  pressure-saturation  relation.  We  also  identify  a  non-toxic  DNAPL 
surrogate  that  can  be  used  in  laboratory  experiments. 

Section  HI  presents  a  unique  method  for  modeling  one-dimensional,  multiphase  flow 
through  layered  media.  This  is  an  important  problem  in  contaminant  hydrology  because 
of  the  prevalence  of  natural  layered  soil  formations  and  the  sensitivity  of  NAPL  flow  to 
small-scale  heterogeneities.  While  it  is  recognized  that  NAPL  transport  is  highly 
sensitive  to  heterogeneities,  currently  available  models  are  poorly  equipped  to 
accommodate  their  influence  on  multiphase  flow.  To  avoid  the  difficulties  found  when 
using  most  other  models,  our  approach  is  to  reduce  the  governing  nonlinear,  partial 
differential  equation  in  time  and  space  to  a  second-order,  ordinary  differential  equation 
by  applying  finite  difference  techniques  to  the  time  derivative  only.  While  this  technique 
has  limited  generality,  it  proved  to  have  significant  advantages  when  modeling  one¬ 
dimensional,  two-phase  flow  through  layered  soils.  By  means  of  this  model,  we  can  gain 
insight  into  multiphase,  flow  processes  at  a  textural  interface. 


In  Section  IV,  we  present  details  of  a  two-dimensional,  pore-scale  model  based  on  a 
modified  diffusion  limited  aggregation  algorithm.  One  of  the  most  stringent  limitations 
of  commonly  used  percolation  models  is  that  they  only  apply  to  low  capillary  number 
flows.  Whereas  this  restriction  is  not  significant  for  LNAPLs,  we  felt  it  was  too  limiting 
for  most  DNAPLs  of  interest  in  contaminant  hydrology.  DNAPLs  are  often  less  viscous 
than  water  and  exhibit  unstable  flow  characteristics.  The  stochastic  aggregation  model 
developed  in  this  project  allows  us  to  more  realistically  characterize  these  unstable  flow 
phenomena. 


The  results  of  this  project  will  be  disseminated  in  refereed  journal  publications.  To  date, 
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SECTION  I 


SECTION  I 


1.1  INTRODUCTION 

Understanding  immiscible  fluid  behavior  in  porous  media,  in  an  effort  to 
characterize  and  remediate  subsurface  sites  contaminated  with  nonaqueous-phase  liquids 
(NAPLs),  has  been  the  focus  of  considerable  research  over  the  last  decade.  This  has  been 
approached  from  many  perspectives  -  from  field  scale  behavior  to  laboratory  studies  to 
analytical  and  numerical  modeling.  Although  research  has  now  generally  progressed  to 
investigation  of  remediation  technologies  and  evaluation  at  the  field  scale,  much  has  been 
learned,  and  remains  to  be  fully  appreciated,  about  the  subtle  aspects  of  immiscible  fluid 
behavior.  A  better  understanding  of  the  small  scale  intricacies  of  NAPL  behavior  can,  and 
has,  led  to  more  accurate  conceptual  models  of  the  subsurface  distribution  and  behavior  of 
these  contaminants.  Our  research  combines  experimental  laboratory  investigations  and 
mathematical  modeling  to  gain  insight  into  the  physical  phenomena  governing  the  flow  and 
distribution  of  NAPLs  in  porous  media.  Emphasis  is  placed  on  the  effects  of  hysteresis  in 
multiphase  flow  in  layered  soils. 

This  report  summarizes  three  different  topics  in  immiscible  fluids.  Section  II 
provides  examples  of  seemingly  subtle  aspects  of  immiscible  fluid  behavior.  However,  it  is 
shown  that  this  behavior  can  be  manifested  as  important  large  scale  effects  in  certain 
scenarios.  Results  of  one  and  two-dimensional  multiphase  laboratory  experiments  of 
transient  flow  in  layered  soils  are  presented.  This  is  an  important  problem  in  contaminant 
hydrology  both  because  of  the  prevalence  of  natural  layered  soil  formations  and  because 
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layers  are  artificially  created  by  processes  such  as  hydrofracturing,  a  practice  that  is 
increasingly  being  used  in  remediation  of  tight  soils.  Capillary  barriers  are  another 
example  where  accurate  simulation  of  multiphase,  transient  flow  across  textural  interfaces 
is  important.  At  a  smaller  scale,  it  is  generally  recognized  that  NAPL  transport  is  highly 
sensitive  to  subsurface  heterogeneities  (Kueper  and  Frind,  1991;  Pantazidou  and  Sitar, 
1993;  Illangasekare,  et  al.,  1995).  Chapter  1  in  section  II  summarizes  a  laboratory 
experiment  conducted  in  a  vertical  column  with  fine  and  coarse  layering.  Water  was 
drained  and  replaced  with  either  oil  or  air.  A  dual  energy  gamma  system  was  used  to 
quantify  fluid  saturations  during  transient  drainage.  The  model  presented  later  in  Section 
III  was  used  to  simulate  this  experiment.  Several  two-dimensional  flume  experiments  are 
also  presented  in  this  section.  The  objective  of  these  experiments  was  to  determine  the 
lens  configuration  of  an  LNAPL  or  DNAPL  lens  in  equilibrium.  Results  show  the 
importance  of  hysteresis  in  the  capillary  pressure-saturation  relation.  We  also  identify  a 
nontoxic  DNAPL  surrogate  that  can  be  used  in  laboratory  experiments. 

Section  III  presents  a  unique  method  for  modeling  one-dimensional,  multiphase 
flow  through  layered  media.  Models  that  simulate  flow  of  immiscible  organic 
contaminants  in  groundwater  are  reviewed  by  Abriola  (1988),  Pankow  and  Cherry  (1996) 
and  Mercer  et  al.  (1996).  However,  as  Mercer  et  al.  note  “immiscible  transport  of 
NAPLs  is  highly  sensitive  to  subsurface  heterogeneities ....  and  the  commonly  available 
models  are  not  well-equipped  to  accommodate  their  influence  in  multiphase  flow”.  To 
avoid  the  difficulties  inherent  in  most  available  models,  our  approach  is  to  reduce  the 
governing  nonlinear,  partial  differential  equation  in  time  and  space  to  a  second-order, 
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ordinary  differential  equation  by  applying  a  finite  difference  approximation  to  the  time 
derivative  only.  Then,  for  each  time  step,  the  resulting  equation  is  solved  as  an  initial 
boundary  problem  by  the  Runge-Kutta  method.  Using  specified  or  assumed  boundary 
conditions,  the  shooting  method  is  used  in  an  iterative  manner  to  meet  the  specific 
boundary  conditions  at  the  other  end  of  the  flow  regime.  This  approach  differs  from 
traditional  finite  difference  or  finite  element  models  in  the  following  ways:  (1)  a  set  of 
equations  is  solved  explicitly  and  sequentially  instead  of  simultaneously,  (2)  the  spacing 
between  nodes  is  varied  within  and  between  time  steps  based  on  accuracy  criteria,  and  (3) 
the  iteration  process  requires  convergence  of  only  one  boundary  condition  value,  rather 
than  iterations  being  performed  simultaneously  on  values  associated  wit  all  internal  nodes. 
The  second  chapter  in  this  section  discusses  the  need  for  multiphase  modeling,  rather  than 
Richards’  equation,  for  simulations  such  as  one-dimensional  drainage,  and  the  third 
chapter  verifies  the  model  by  comparing  simulations  with  laboratory  data. 

Section  IV  presents  details  of  a  two-dimensional  pore  scale  model  based  on  a 
modified  diffusion  limited  aggregation  algorithm.  We  also  developed,  as  part  of  this 
project,  an  invasion  percolation  model  to  simulate  pore  scale  processes  in  porous  media, 
and  predict  macroscopic  capillary  pressure-saturation  relations  and  relative  hydraulic 
conductivity  for  strongly  water-wetted  systems  at  low  capillary  numbers.  After 
completion  of  the  percolation  model,  however,  we  revised  our  approach  to  pore  scale 
modeling.  One  of  the  most  stringent  limitation  on  models  based  on  percolation  is  that  they 
are  limited  to  low  capillary  numbers.  Whereas  this  restriction  is  not  significant  for 
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LNAPLs,  we  felt  it  was  too  limiting  for  modeling  dense  non-aqueous  phase  liquids  which 
exhibit  unstable  flow. 
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SECTION  II 


SECTION  II 


Chapter  1 

Experimental  Study  of  One-Dimensional  Immiscible  Fluid  Drainage 

In  Layered  Sands 


Abstract 

Data  that  characterize  transient  flow  of  nonaqueous  phase  liquids  in 
heterogeneous  media  are  needed  for  verification  of  multiphase  flow  models.  This  paper 
summarizes  three  drainage  experiments  in  layered  sands  performed  specifically  to 
provide  such  data.  The  layered  system  consists  of  a  coarse  sand  between  two  finer  sands. 
In  all  cases,  water  was  the  draining  fluid.  In  the  first  experiment,  air  at  atmospheric 
pressure  was  allowed  to  move  freely  into  and  out  of  the  coarse  sand  layer.  In  the  second 
experiment,  air  entry  into  the  coarse  layer  was  only  possible  through  the  overlying  fine 
sand.  These  two  experiments  illustrated  effects  of  nonwetting  phase  accessibility.  In  the 
third  experiment,  water  was  displaced  by  a  light  nonaqueous  phase  liquid.  Both  the  effect 
of  nonwetting  phase  access  to  the  coarse  material  and  the  effect  of  a  non-negligible 
viscosity  were  evident  in  this  experiment.  The  three  experiments  represent  a  typical 
progression  in  multiphase  model  development.  The  layered  system  results  in  flow 
conditions  that  give  many  conventional  models  computational  problems  due  to  property 
discontinuities  at  layer  interfaces. 


1.1  INTRODUCTION 


The  scale  at  which  subsurface  heterogeneities  can  be  characterized  in  a  field 
situation  is  limited.  However,  a  better  understanding  of  how  heterogeneities  affect  the 
transport  and  distribution  of  nonaqueous  phase  liquids  (NAPLs)  on  a  small  scale  leads  to 
more  accurate  conceptualizations  of  site  contamination,  better  interpretations  of  observed 
field  data,  and  more  realistic  models.  The  conceptual  model  assumed  for  the  distribution 
of  a  contaminant  greatly  affects  the  determination  of  important  site  characteristics  such  as 
rate  of  mass  transfer  by  diffusion,  mobility  of  a  contaminant,  the  areal  extent  of 
contamination,  and  even  the  estimation  of  total  mass  present.  Therefore,  understanding 
small  scale  intricacies  of  NAPL  transport  is  relevant  to  interpretation  of  field  data  and 
design  of  field  scale  remediation  strategies. 

In  modeling  transient,  multiphase  flow  with  conventional  finite  difference  and 
finite  element  models,  numerical  instabilities  and  mass  balance  errors  can  occur  if  there 
are  saturation  discontinuities  and/or  large  pressure  gradients  (Celia  et  al.,  1990;  Segol, 
1994).  The  problem  is  particularly  acute  when  a  sharp  saturation  front  crosses  the 
interface  between  two  dissimilar  soils.  Nonetheless,  the  importance  of  adequately 
characterizing  layered  or  otherwise  heterogeneous  media  has  been  clearly  shown  in  many 
experimental  studies  (e.g.,  Kueper  and  Frind,  1991;  Pantazidou  and  Sitar,  1993; 
Illangasekare  et  al.,  1995a,  1995b).  The  objective  of  this  work  was  to  provide  simple 
laboratory  data  for  model  verification  of  two-phase  flow  in  layered  porous  media. 
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1.2  MATERIALS  AND  METHODS 


Experiments  were  conducted  in  a  vertical  column  packed  with  fine  sand  from  the 
bottom  of  the  column  to  a  height  of  43  cm,  a  coarse  sand  from  43  to  57.5  cm,  and  fine 
sand  again  from  57.5  cm  to  80  cm.  Three  experiments  were  conducted.  In  experiment 
one,  water  drained  through  the  bottom  of  the  column  and  was  replaced  with  air  that  Was 
allowed  to  freely  enter  and  exit  the  column  through  ports  at  the  top  of  the  coarse  sand 
layer.  Valves  allowing  air  entry  were  opened  as  soon  as  the  water  pressure  in  the  coarse 
layer  became  less  than  atmospheric.  In  the  second  experiment,  the  air  was  again  the 
displacing  fluid  but  free  access  was  not  provided  into  and  out  of  the  column  through  the 
ports.  The  only  pathway  for  the  air  was  through  the  overlying  finer  layer.  In  the  third 
experiment,  a  lighter-than-water  nonaqueous  phase  liquid  (LNAPL)  was  the  displacing 
fluid. 

The  inside  dimensions  of  the  experimental  column  were  91  cm  high,  15.3  cm 
wide  and  5.6  cm  thick.  The  front  and  back  walls  were  constructed  with  0.635  cm-thick 
glass.  The  side  walls  were  aluminum.  Twelve  pressure  ports  were  drilled  into  one  side 
of  the  column,  two  each  at  six  locations  along  the  height  of  the  column.  Hydrophilic  and 
hydrophobic  pressure  tensiometers  (Noller,  1992)  were  installed  at  each  of  the  six 
locations  and  connected  to  a  data  acquisition  system. 

Fluid  saturations  were  measured  using  a  dual-energy  gamma  attenuation  system. 
Gamma  attenuation  theory  and  its  application  to  multiphase  flow  problems  in  porous 
media  are  well  documented  in  the  literature,  including  discussions  of  the  errors  typically 
observed  with  gamma  attenuation  measurements  (Eckberg  and  Sunada,  1984;  Schiegg 
and  McBride,  1987;  Nofziger  and  Swartzendruber,  1974;  Stroosnijder  and  DeSwart, 
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1 974;  Lenhard  et  al.,  1988,  1991 ;  Lenhard  and  Parker,  1988;Ferrandetal.,  1986,  1989; 

Illangasekare  et  al.,  1995b;  Oostrom  et  al.,  1995).  Details  of  the  system  used  in  this  study 
are  given  by  Noller  (1992). 

The  porous  media  used  for  this  investigation  were  two  different  clean  sands 

obtained  from  Colorado  Silica  Sand,  Inc.,  Colorado  Springs,  Colorado.  The  finer  sand  is 

referred  to  as  Colorado  Springs  Silica  Sand  (CSSS)  #100  and  the  coamer  sand  is  referred 

to  as  CSSS  #40-60.  Physical  and  hydraulic  properties  of  the  sands  are  shown  in  Table 

1 . 1 .  Both  Brooks-Corey  (Brooks  and  Corey,  1 966)  and  van  Genuchten  (van  Genuchten, 

1980)  parameters  were  fit  to  water-air  retention  data  for  the  two  sands  and  are  also  given 

tn  Table  1.1.  (The  water  retention  curves  and  particle  size  distributions  are  included  in 

Appendix  A.)  The  column  was  packed  following  the  mining  method  of  Rad  and  Tumay 
(1985). 


Table  1.1.  Sand  properties. 


Dpo  (mm) 

Uniformity  Coefficient  (d^/d,,,) 
Brooks-Corey  P j  (cm  of  water) 
Brooks-Corey  X 
van  Genuchten  n 
van  Genuchten  ot  (cm  of  water)-1 
Residual  Water  Saturation,  Sw 
Avg.  Packed  Porosity 
Hydraulic  Conductivity  (cm  s  ') 


CSSS  #100 
0.16 
0.08 
0.30 
2.18 
41 
1.7 
4.4 
0.016 
0.08 
0.37 
0.002 


0.25 

0.30 

1.12 

18 

3.4 

6.3 

0.041 

0.05 

0.37 

0.010 
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Distilled  water  was  used  as  the  wetting  fluid  in  the  column  experiments.  The 
nonwetting  fluids  were  air  and  a  15:1  ratio  by  volume  of  Soltrol  220  and  1-iododecane, 
CH3(CH2)9 1.  Soltrol  220  is  a  clear,  isoparaffinic  solvent  manufactured  by  Phillips 
Petroleum  Co.,  Bartlesville,  Oklahoma.  The  presence  of  iodide  increases  the  gamma 
attenuation  by  the  Soltrol  and  enhances  differences  in  the  attenuation  coefficients  of  the 
water  and  Soltrol  (Lenhard  et  al.,  1988).  An  organic  dye,  Sudan  IV  (Aldrich  Chemical 
Co.,  Milwaukee,  Wisconsin)  was  added  to  the  Soltrol  to  allow  visual  differentiation 
between  the  LNAPL  and  the  water  phases.  Properties  of  the  Soltrol/iododecane/dye 
mixture  are  given  in  Table  1 .2.  Before  each  experiment,  the  soil  column  was  flushed 
with  C02  and  then  saturated  with  distilled  water  from  the  bottom  at  a  rate  of  3  ml/min. 
An  initial  wetting  phase  saturation  of  approximately  90%  was  obtained  for  all 
experiments. 


Table  1.2.  Fluid  properties  of  water  and  Soltrol/iododecane/dye  mixture  (at  22°C). 


water 

LNAPL 

Density  (g  cm-3) 

0.997 

0.83 

Viscosity  (centipoise) 

0.95 

3.65 

Surface  Tension  (dynes  cnrl) 

70 

28 

Interfacial  Tension  (dynes  cm-1) 

N/A 

42a 

Gamma  Attenuation  Coeff.  (cm2  g-1) 

24lAm 

0.205 

0.508 

l37Cs 

0.092 

0.088 

*  values  decrease  from  42  to  35  dynes  cm'1  with  water  contact.  The  model  presented  later  used  the  initial  value  42  dynes  cm'1. 
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1.3  experimental  procedures 

Experiment  1:  Water  drainage  with  “inviscid  air- at  atmospheric. pressure.  A 
firs,  step  in  multiphase  model  verification  is  .omodei  an  air-water  flow  problem 
assuming  the  air  can  move  freely  into  and  out  of  the  flow  domain  with  the  pressure  of  the 
air  everywhere  atmospheric.  In  experiment  1 ,  this  scenario  was  simulated  by  providing 

inlet  into  the  coarse  layer.  This  was  recently  described  by  Marinelli  (1996)  as  the 
more  accurate  view  of  the  air  phase  in  the  Richards  efiuation  not  same  as 

assuming  the  air  to  be  inviscid,  i,„  having  zero  viscosity.  The  column  was  initially 

distilled  water.  At  time  t  =  0,  the  water  table  was  dropped  instantaneously 
from  just  above  the  top  of  the  soil  to  an  elevation  10  cm  below  the  bottom  of  the  column 

Thus,  tile  boundary  conditions  for  experiment  1  are  a  water  pressure  of -10  cm  of  water  at 

the  bottom  of  tire  column  and  atmospheric  air  pressure  at  the  top  of  the  column.  Data 
were  collected  over  57  hours. 

Experiment  2:  Water  drainage  with  nir  r. 

mage  with  air.  In  expenment  2,  the  initial  and 

boundary  conditions  of  experiment  1  were  repeated.  However  in  this  • 

nowever,  m  this  experiment,  the  air 

was  not  provided  a  free  pathway  inti)  the  middle  iayer.  Air  couid  only  access  the  middle 
coame  layer  through  the  overlying  finer  soil.  This  may  be  mom  realistic  in  a  field 
situation  and  represents  a  progression  in  modeling  difficulty.  I„  this  case,  water  a, 
pressures  less  than  atmospheric  in  die  middle  coarse  layer  will  not  drain  until  a  pathway 
through  the  finer  overlying  soil  exists.  Dam  were  collected  over  34  hours 
Crimen,  3:  Water  drainage  Ki,h  an  LNAPL.  experiment  3,  an  LNAPL 
was  the  displacing  fluid.  Again,  the  column  was  initially  sa, mated  with  water,  with  a 
small  ponded  depth  of  water  just  above  the  top  of  the  soil  column.  Approximately  5  cm 


of  LNAPL  was  added  on  top  of  the  water.  At  time  t  =  0,  the  water  table  was  dropped 
instantaneously  to  an  elevation  10  cm  below  the  bottom  of  the  column.  The  oil  table  was 
maintained  at  5  cm  above  the  top  of  the  column.  Therefore,  the  boundary  conditions  for 
this  experiment  were  a  water  pressure  of  -10  cm  of  water  at  the  bottom  of  the  column, 
and  an  oil  pressure  of  5  cm  of  oil  at  the  top  of  the  column. 

Drainage  was  allowed  to  continue  until  the  oil-water  interface  approached  the 
bottom  of  the  column,  and  was  stopped  by  changing  the  lower  boundary  pressure  to  80 
cm  of  water  at  time  t  =  375  minutes.  The  subsequent  reimbibition  of  water  was 
monitored  for  27  hours. 


1.4  RESULTS  AND  DISCUSSION 

Fluid  saturations  and  outflow  rates  of  the  wetting  fluid  were  recorded  and  yielded 
good  information.  Pressure  data  was  recorded  for  the  water-air  experiments  but,  due  to 
equipment  failures,  only  limited  pressure  data  was  taken  for  the  water-oil  experiment. 
Summaries  of  these  data  are  presented  here.  The  set  of  outflow  and  saturation  data  is 
included  in  Appendix  A. 

Figure  1.1  shows  wetting  phase  outflow  rate  versus  time  for  the  three 
experiments.  In  both  water-air  drainage  experiments,  the  coarse  sand  layer  desaturated 
while  saturations  in  the  upper  sand  layer  were  still  quite  high.  The  effect  of  this  early 
desaturation  of  the  coarse  layer  is  shown  by  an  abrupt  decrease  in  outflow  rate  in 
experiments  1  and  2  between  50  and  80  minutes.  When  the  coarse  layer  desaturated,  the 
relative  permeability  in  that  layer  decreased  significantly,  retarding  drainage  of  the  water 
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from  the  overlying  finer  material.  The  low  relative  permeability  in  the  coarse  material 

and  consequent  reduction  in  drainage  rate  illustrate  the  mechanism  behind  capillary 
barriers. 


Figure  1.1  Outflow  rate  with  time  for  experiments  1, 2,  and  3. 


The  water  pressure  in  the  coarse  layer  can  become  less  than  atmospheric  without 
desaturation  of  the  pores  if  air  entry  is  not  permitted.  If  the  middle  layer  stays  saturated 
longer,  the  equivalent  overall  permeability  will  be  larger  for  a  longer  time  and  a 
difference  should  be  seen  in  the  column  outflow  rate.  Differences  in  the  drainage  rate 
between  experiments  1  and  2  in  figure  1.1  are  not  evident  when  plotted  on  the  same  curve 
as  the  water-LNAPL  experiment  data.  This  is  because  the  viscosity  of  the  air  is  about  50 
times  less  than  that  of  water,  so  when  there  is  even  a  slight  desaturation  of  the  fine  layer, 
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the  air  flows  freely  through  the  surface  layer  and  the  middle  layer  desaturates  quickly. 

The  LNAPL  has  a  viscosity  four  times  larger  than  water  so  even  when  access  is  available 
for  the  nonwetting  fluid,  the  coarse  layer  doesn’t  drain  as  quickly.  Therefore,  differences 
between  the  air  and  LNAPL  nonwetting  phases  are  much  more  evident  in  figure  1.1  than 
differences  between  the  two  air  experiments.  However,  differences  between  experiments 
1  and  2  can  be  seen  by  comparing  saturation  profiles  and  cumulative  outflow  from  these 
two  experiments.  Figure  1.2  shows  wetting  phase  saturations  at  the  same  time  from  the 
two  water-air  experiments,  with  higher  saturations  remaining  in  experiment  1.  Figure  1.3 
shows  that  cumulative  outflow  in  experiment  1  lags  behind  cumulative  outflow  in 
experiment  2.  As  long  as  the  coarse  layer  stays  saturated,  the  wetting  phase  relative 
permeability  in  the  coarse  layer  equals  the  saturated  permeability  for  that  layer.  Drainage 
is  slowed  when  the  middle  layer  desaturates  and  permeability  decreases.  The  cumulative 
outflow  graph  also  clearly  shows  a  change  in  slope  for  each  of  the  water-air  experiments. 
These  slopes  correspond  to  the  drainage  rate  before  and  after  the  coarse  layer  desaturated. 


Figure  1.2.  Volumetric  water  content  with  elevation  for  air  experiments  at  3.5  hours. 
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The  water-LNAPL  experiment  shows  effects  both  of  accessibility  and  of  a 
nonnegligible  viscosity  of  the  nonwetting  phase.  In  experiment  three,  the  LNAPL  is  not 
provided  free  access  to  the  middle  coarse  layer  and,*in  addition,  the  LNAPL  viscosity  is 
about  four  times  larger  than  the  water  viscosity.  The  result  is  the  water-LNAPL  front 
moves  from  the  top  of  the  column  downward  at  a  relatively  uniform  rate.  This  is  shown 
in  figure  1.4.  The  middle  layer  cannot  desaturate  until  the  LNAPL  has  both  access  to  the 
layer  through  the  overlying  fine  soil  and  time  to  move  into  the  coarser  layer.  The  rate 
that  the  LNAPL  can  move  through  the  fine  layer  to  the  coarse  layer  is  dependent  on  the 
NAPL  saturation  in  the  fine  layer.  If  the  NAPL  saturation  is  low,  i.e.,  the  fine  layer  has 
not  drained  significantly,  the  permeability  of  the  NAPL  in  the  top  layer  is  small. 
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Figure  1.4.  Migration  of  the  LNAPL  front. 


There  are  some  subtle  features  of  the  figures  presented  that  are  worth  noting  and 
lend  credibility  to  the  gamma  data.  In  figure  1.4,  the  highest  saturation  of  the  LNAPL  in 
the  coarse  layer  is  slightly  higher  than  the  saturation  of  the  LNAPL  in  the  overlying  fine 
layer.  This  reflects  a  higher  residual  water  content  in  the  fine  layer.  In  figure  1 .3,  there 
are  slight  changes  in  slope  noticeable  if  early  data  and  late  data  are  compared  for 
experiment  3.  The  higher  slope  on  the  cumulative  outflow  curve  would  be  expected 
because  relatively  more  of  the  column  being  drained  at  this  time  is  coarse  material  with 
higher  saturated  permeability. 
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Another  interesting  feature  is  an  abrupt  decrease  in  outflow  rate  in  experiment  3  at 
about  200  minutes.  This  is  shown  clearly  in  figure  1.5.  From  the  saturation  data  (figure 
1 .4),  it  can  be  seen  that  the  abrupt  drop  in  outflow  rate  corresponds  to  the  time  when  the 
water-LNAPL  interface  moved  from  the  coarse  layer  into  the  lower  fine  layer.  This 
effect  is  also  captured,  at  least  qualitatively,  by  the  numerical  modeling  of  this 
experiment  also  shown  in  figure  1.5.  It  should  be  noted  that  the  model  shown  in  figure 
1.5  used  the  parameters  presented  in  tables  1.1  and  1.2  with  no  calibration  for  a  fit.  For 
more  information  on  the  numerical  model  and  a  similar  experiment  the  reader  is  referred 
to  Marinelli  (1996).  The  capillary  pressure  required  for  pore  drainage  is  inversely  related 
to  pore  size  by  the  LaPlace  equation  of  capillarity  (Corey,  1994).  When  the  interface 
moved  into  the  fine  layer,  the  pore  size  decreased  and  the  capillary  pressure  had  to 
increase  before  drainage  could  continue.  Thus,  the  gradient  also  decreased,  reducing 
outflow.  Since  the  fluid-fluid  front  is  sharp,  the  drop  in  outflow  is  also  abrupt. 

At  time  t  =  375  minutes,  drainage  in  experiment  3  was  stopped  by  changing  the 
lower  boundary  condition  to  a  positive  water  pressure  of  80  cm,  i.e.,  raising  the  water 
table  to  the  top  of  the  column.  This  reversed  the  gradient  and  water  reimbibed  into  the 
column,  approaching  an  equilibrium  distribution  over  the  next  24  hours.  Figure  1.6 
shows  saturations  at  time  t  =  33  hours.  This  figure  shows  high  LNAPL  contents  just 
below  the  interface  between  the  coarse  and  fine  layers.  This  trapping  of  LNAPL  below  a 
finer  layer  has  been  observed  in  field  studies  (Vroblesky  et  al.,  1995). 
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Figure  1.5.  Outflow  rate  for  the  LNAPL  experiment  showing  an  abrupt  decrease  at  200  minutes. 


Figure  1.6.  Fluid  contents  at  33  hours,  27  hours  after  the  water  table  was  returned  to  the  top  of  the 
column. 
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1.5  CONCLUSIONS 


The  laboratory  experiments  presented  in  this  paper  provide  data  for  drainage  of 
immiscible  fluids  in  layered  soils.  Soil  layering  results  in  discontinuities  in  saturations 
and  permeabilities,  effects  that  often  create  numerical  difficulties  in  conventional  finite 
difference  and  finite  element  models.  The  data  acquired  are  sufficiently  detailed  to  be 
useful  for  evaluation  of  two-phase  flow  models.  Effects  of  soil  layering,  nonwetting 
phase  access,  and  differences  in  those  effects  for  air  and  LNAPL  systems  were  observed. 


1.6  REFERENCES 

Brooks,  R.  H.  and  A.  T.  Corey.  1966.  Properties  of  porous  media  affecting  flow. 

Journal  of  Irrigation  and  Drainage  Div.,  ASCE,  92(IR2):61-68. 

Celia,  A.  M.,  E.T.  Bouioutas  and  R.L.  Zarba.  1990.  A  general  mass-conservation 

26m"?83-°4960n  UnSa,Urated  flow  e9“a«on.  Water  Resources  Research 

Corey,  A.  T.  1994.  Mechanics  of  Immiscible  Fluids  in  Porous  Media.  Water 

Resources  Publications,  Littleton,  CO,  252  pp. 

Eckberg,  D.K.  and  D.K.  Sunada.  1984.  Nonsteady  three-phase  immiscible  fluid 

is  n  ution  in  porous  media.  W, ater  Resources  Research.  20: 1 89 1  - 1 897 

FenandLARC.D.  Milly,  and  G.T.  Pinder.  1986.  Dual-gamma  attenuation  for  the 

Wi“' KSPeC' ,0  ^  flUidS' 

Ferrand,DA.,  P.CD.  Milly,  and  O.F.  Finder.  1989.  Experimental  determination  of 
*££%££**  “  P°r0US  media'  J°mnal  <>f  Contaminant 

Illangasekare,  T.H.,  E.J.  Atmbruster  III  and  D.N.  Yates.  1995a.  Nonaqueous-phase 

**■  J°“™'  of  Environmental 


11-14 


Illangasekare,  T.H.,  J.L.  Ramsey,  Jr.,  K.H.  Jensen,  M.B.  Butts.  1995b.  Experimental 
study  of  movement  and  distribution  of  dense  organic  contaminants  in 
heterogeneous  aquifers.  Journal  of  Contaminant  Hydrology.  20:1-25. 

Kueper,  B.  H.  and  E.O.  Frind.  1991.  Two-phase  flow  in  heterogeneous  media.  1 .  Model 
development.  Water  Resources  Research.  27(6):1049-1057. 

Lenhard,  R.J.,  J.H.  Dane,  J.C.  Parker,  and  J.J.  Kaluarachchi.  1988.  Measurement  and 
simulation  of  one-dimensional  transient  three-phase  flow  for  monotonic  liquid 
drainage.  Water  Resources  Research.  24:853-863. 

Lenhard,  R.J.  and  J.C.  Parker.  1988.  Experimental  validation  of  the  theory  of  extending 
two-phase  saturation-pressure  relations  to  three-fluid  phase  systems  for  monotonic 
drainage  paths.  Water  Resources  Research.  24:373-380. 

Lenhard,  R.J.,  J.C.  Parker  and  J.J.  Kaluarachchi.  1991 .  Comparing  simulated  and 
experimental  hysteretic  two-phase  transient  fluid  flow  phenomena.  Water 
Resources  Research.  27:2113-2124. 

Marinelli,  F.  1996.  Semi-analytical  solutions  for  one-dimensional  flow  of  immiscible 
fluids  in  layered  porous  media.  Ph.D.  dissertation.  Colorado  State  University, 

Fort  Collins,  CO.  124  pp. 

Nofziger,  D.L.  and  D.  Swartzendruber.  1974.  Material  content  of  binary  physical 
mixtures  as  measured  with  a  dual-energy  beam  of  gamma  rays.  Journal  of 
Applied  Physics.  45:5443-5449. 

Noller,  K.L.  1992.  Calibration  of  a  dual-beam  gamma  system  to  investigate  multiphase 
equilibrium  profiles  in  soil.  Master’s  thesis.  Colorado  State  University,  Fort 
Collins,  CO.  163  pp. 

Oostrom,  M.,  J.H.  Dane,  B.C.  Missildine  and  R.J.  Lenhard.  1995.  Error  analysis  of  dual¬ 
energy  gamma  radiation  measurements.  Soil  Science.  160(l):28-42. 

Pantazidou,  M.,  and  N.  Sitar.  1993.  Emplacement  of  nonaqueous-phase  liquids  in  the 
vadose  zone.  Water  Resources  Research.  29(3):705-722. 

Rad,  N.S.  and  M.  T.  Tumay.  1985.  Factors  affecting  sand  specimen  preparation  by 
raining.  Geotechnical  Testing  Journal.  ASTM.  10(1):3 1  -37. 

Scheigg,  H.O.  and  J.F.  McBride.  1987.  Laboratory  setup  to  study  two-dimensional 

multiphase  flow  in  porous  media.  Proceedings  of  the  NWWA/API  Conference  on 
Petroleum  Hydrocarbons  and  Organic  Chemicals  in  Ground  Water  -  Prevention, 
Detection  and  Restoration,  Houston,  Texas,  pp.  371-394. 


11-15 


’  nc->  Englewood  Cliffs, 

ktroosnijder,  L.  and  J  G  rvc 

Vroblesky,  D  A  J  F  R  5 


11-16 


Chapter  2 

Geometry  of  NAPL  Mounds  in  Water-NAPL  Systems 


Abstract 

Nonaqueous-phase  liquids  (NAPLs)  in  the  subsurface  can  be  distributed  in 
mounds  or  pools  of  significant  size  and  thickness,  at  relatively  high  NAPL  saturations.  A 
depression  in  a  fine  textured  heterogeneity  is  not  required  for  this  to  occur.  The 
mechanics  of  immiscible  fluids  allow  lateral  stability  even  on  perfectly  flat  or  uniformly 
sloped  strata.  This  paper  summarizes  laboratory  experiments  performed  to  illustrate  this 
behavior  and  to  verify  an  expression  given  to  estimate  the  maximum  height  of  a  static 
NAPL  mound.  It  is  shown  that  this  analysis  applies  both  to  DNAPLs  (denser-than-water 
NAPLs)  on  top  of  heterogeneities  and  to  air  and  LNAPLs  (lighter-than- water  NAPLs) 
trapped  below  the  water  table  beneath  heterogeneities.  It  is  also  shown  that  this  behavior 
cannot  be  explained  without  considering  the  hysteresis  in  capillary  pressure-saturation 
relationships.  Models  not  incorporating  this  hysteresis  may  result  in  incorrect 
simulations  of  the  migration  and  distribution  of  NAPLs  in  the  subsurface. 

2.1  INTRODUCTION 

The  behavior  of  NAPLs  in  the  subsurface  has  been  the  subject  of  much  research 
in  recent  years.  Groundwater  contamination  from  NAPLs  such  as  gasoline  and 
chlorinated  solvents  is  a  widespread  problem.  An  understanding  of  the  principles 


11-17 


involved  with  NAPLs  in  groundwater  is  necessaiy  for  accurate 


assessment  and  effective 


remediation  of  these  sites. 


The  transport  and  distribution  of  NAPLs  is  strongly  affected  by  subsurface 

heterogeneities.  Although  the  detail  to  which  subsurface  hetetogeneities  can  be 

characterized  in  the  field  is  limited,  a  better  understanding  of  how  these  heterogeneities 

affect  NAPL  distribution  should  lead  to  more  accurate  conceptualizations  of  site 

contamination,  better  interpretations  of  observed  field  data,  and  more  realistic  models. 

The  conceptual  model  assumed  for  the  distribution  of  a  free-phase  contaminant  greatly 

affects  the  estimation  of  important  site  characteristics  such  as  the  rate  of  mass  transfer  by 

diffusion  or  dissolution,  the  mobility  of  fire  contaminant,  and  fire  areal  extent  of  a 
contaminant. 


Pools  or  mounds  of  denser-than-water  NAPLs  (DNAPLs)  commonly  occur  on 
low  permeability  lenses  or  strata  in  fire  subsurface.  The  trapping  of  lighter-flrarr-water 
NAPLs  (LNAPLs)  under  lenses  below  fire  water  table  has  also  been  observed  wifi,  jet 

“  (Vr0WeSky  e‘ al  ’  1995)  “■  *“■  *  *  sParfihig  (Lundegard  and  Anderson, 

1996).  Depressions  in  low  permeability  lenses  are  no,  teqnired  for  NAPL  pools  of 

significant  size  to  form,  they  can  torn,  on  flat  lenses  and  even  in  homogeneous  material 
(Keuper  e,  al.,  1993).  Understanding  the  behavior  and  distribution  of  these  pools 
contributes  to  the  characterization  and  remediation  of  the  numerous  DNAPL 
contaminated  sites. 


McWhorter  and  Kueper  (1996)  and  Kueper  et  al.  (1993)  explain  how  the 
mechanics  involved  with  NAPL  pools  can  be  used  to  interpret  observed  DNAPL 
thicknesses  in  monitoring  wells,  and  to  estimate  the  thickness  and  volume  of  DNAPL 
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pools.  Lundegard  and  Anderson  (1996)  discuss  air  mounding  beneath  lenses  during  air 
sparging.  They  mention  this  mounding  can  lead  to  complex  air  flow  patterns  and  can 
greatly  increase  lateral  groundwater  displacement.  This  exaggerated  displacement  under 
heterogeneous  conditions  casts  doubt  on  its  use  as  an  indicator  of  the  air  sparging  radius- 
of-influence.  They  conclude  that  under  heterogeneous  conditions,  detailed  knowledge  of 
site  stratigraphy  and  sediment  properties  is  necessary  to  predict  the  flow  response  of  air 
injection.  Accurately  estimating  the  extent  of  air  trapped  under  lenses  could  improve 
these  predictions. 

The  objectives  of  this  study  were  1)  to  predict  the  maximum  thickness  of  laterally 
stable  NAPL  mounds  by  accounting  for  hysteresis  in  the  capillary  pressure-saturation 
(Pc-S)  relationships,  2)  to  extend  this  analysis  to  LNAPL  and  air  mounds  trapped  beneath 
low  permeability  lenses,  and  3)  to  provide  laboratory  verification  of  the  predictive 
equations  of  this  immiscible  fluid  behavior.  It  is  shown  that  the  phenomena  of  static 
NAPL  mounds  on  subsurface  heterogeneities  cannot  be  fully  explained  without 
considering  hysteresis. 

2.2  CAPILLARY  PRESSURE-SATURATION  RELATIONSHIPS 

In  most  aquifer  materials,  NAPLs  are  nonwetting  with  respect  to  water,  meaning 
water  has  the  greater  attraction  to  the  matrix.  Combined  with  a  non-zero  interfacial 
tension  between  two  immiscible  fluids,  this  enables  a  pressure  difference  to  be  supported 
in  porous  media  (Corey,  1994).  This  pressure  difference  is  termed  the  capillary  pressure 
and  is  defined  as  the  nonwetting  fluid  pressure  minus  the  wetting  fluid  pressure.  If  the 
capillary  pressure  is  increased,  the  nonwetting  fluid  (the  NAPL  or  air)  can  displace  water 
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from  the  matrix.  Figure  2.1  shows  a  typical  relationship  between  capillary  pressure  and 
NAPL-water  saturation  for  a  volume  of  porous  media. 

The  nonwetting  phase  entry  pressure  of  a  porous  medium  refers  to  the  capillary 
pressure  below  which  little  or  no  desaturation  of  the  water  occurs  (Corey,  1994).  Above 
this  pressure,  the  water  in  some  pores  begins  to  be  displaced  by  the  nonwetting  phase, 
and  saturation  decreases  rapidly  along  the  primary  drainage  curve,  eventually 
approaching  an  irreducible  or  residual  level  of  water  saturation.  If  the  capillary  pressure 
is  then  decreased,  the  water  will  begin  to  reimbibe  into  the  media.  However,  hysteresis  is 
exhibited  in  the  capillary  pressure-saturation  (Pc-S)  relationship.  As  the  water  reimbibes 
from  a  residual  saturation,  the  saturation  increases  along  a  secondary  curve,  termed  the 
primary  wetting  or  imbibition  curve.  As  the  capillary  pressure  goes  to  zero,  not  all  the 
nonwetting  fluid  is  displaced,  leaving  a  trapped  or  residual  nonwetting  phase. 


Figure  2. 1 .  A  typical  relationship  between  capillary  pressure  and  fluid  saturation  in  a  volume  of  porous 
media,  illustrating  the  hysteresis  exhibited  between  drainage  and  imbibition  curves. 
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Although  much  less  used  and  documented  in  the  literature,  the  imbibition  curve  of 
the  Pc-S  relationship  often  exhibits  an  inflection  point  roughly  analogous  to  entry 
pressure  of  the  drainage  curve  (Baker  and  Hillel,  1990;  Liu  et  al.,  1993;  Hogarth  et  ah, 
1988).  Below  this  point,  water  saturation  changes  only  slightly  with  decreasing  capillary 
pressure.  Particularly  in  a  uniform  sand,  it  is  near  this  pressure  that  NAPL  saturations 
dramatically  decrease  and  water  saturations  increase  with  decreasing  capillary  pressures. 
This  inflection  point  is  referred  to  here  as  the  imbibition  entry  pressure. 


2.3  THEORY  OF  NAPL  MOUND  FORMATION 

DNAPLs  that  are  spilled  or  leaked  into  the  subsurface  move  through  the  vadose 
zone  under  the  influence  of  gravity  and  capillary  forces.  If  sufficient  volume  is  spilled  a 
DNAPL  will  reach  and  penetrate  the  water  table.  Depending  on  the  volume  and  rate  of 
the  spill,  when  a  DNAPL  migrating  downward  below  the  water  table  encounters  a  fine 
lens,  it  can  begin  to  pool,  mound  up,  and  spread  across  the  lens.  Once  the  source  of  the 
DNAPL  spill  is  exhausted,  the  mound  height  will  decrease  as  the  DNAPL  spreads 
laterally.  However,  as  asserted  by  Kueper  et  al.  (1993),  a  mound  will  not  spread 
indefinitely,  even  on  a  flat  lens.  The  lateral  spreading  of  a  NAPL  collecting  and 
mounding  on  even  a  flat  lens  will  stop  when  the  mound  height  has  decreased  to  a  point 
where  there  is  no  longer  sufficient  capillary  pressure  at  the  advancing  edge  to  overcome 
the  entry  pressure  of  the  aquifer  material.  The  maximum  thickness  of  this  NAPL  mound 
is  dependent  on  the  density  difference  of  the  two  fluids  and  the  entry  pressure  of  the 
surrounding  matrix.  These  parameters  can  lead  to  significant  volumes  of  NAPL 
associated  with  heterogeneities. 
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Figure  2.2  shows  the  hydrostatic  pressure  distribution  in  a  DNAPL  mound.  since 
0  DNAPL  IS  denser  than  water,  its  pressure  increases  faster  with  depth  than  the  water 
pressure.  Therefore,  the  cap, Mary  pressute  increases  with  depth  a  maximum  a.  the 

bottom  of  the  tnonnd.  In  a  sitniiar  anaiysis,  Kueper  e,  ah  („93)  and  McWhorter  <nJ 

Kueper  ( 1 996)  explain  that  for  the  mound  to  remain  stationary,  the  capillary  pressure  at 
the  bottom  of  the  mound  has  to  be  equal  to  or  less  than  the  entry  pressure  of  die  aquitard 
matena,  upon  which  it  is  perched.  This  is  used  estimate  the  maximum  poo,  thickness 
before  penetration  into  the  aquitard.  However,  for  the  point  of  incipient  lateral  motion, 
capillary  pressure  at  the  bottom  leading  edge  of  the 


entry  pressure  of  the  surrounding  media. 


1  mound  (Point  1)  is  equal  to  the 


saturated  po^s^nedL^Pe^enotes  water  dra'*”  DNAPL  mound  perehed  on  ,  fine  [extured  |e„s 

pressure  of  the  media.  “*  mKr  *■"”«»  mry  pressure  and  Pe.,  the  w„«r  imbibid„„  ent^ 
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From  this  analysis  the  capillary  pressure  at  all  fluid-fluid  interfaces  must  be  less 
than  the  nonwetting  phase  entry  pressure  of  the  media  except  at  the  bottom  of  the  mound. 
This  would  seem  to  imply  that  reinvasion  of  the  water  would  occur  and  therefore  the 
DNAPL  could  not  exist  in  mounds.  This  is  resolved  by  considering  the  hysteresis  in  the 
Pc-S  relationships.  The  top  of  the  DNAPL  mound  (Point  2  in  Figure  2.2)  would  be 
receding,  i.e.,  DNAPL  draining  from  the  pores.  Therefore,  the  saturations  at  the  top 
would  be  given  by  the  secondary  Pc-S  curve,  the  water  imbibition  curve,  if  we  assume 
the  water  in  the  DNAPL  invaded  zone  had  decreased  close  to  residual  contents.  This 
would  allow  high  DNAPL  saturations  to  remain  at  low  capillary  pressures.  The  upper 
boundary  of  the  mound,  where  DNAPL  contents  rapidly  decrease  from  high  to  residual 
saturations,  would  correspond  to  the  point  where  the  capillary  pressure  is  equal  to  the 
imbibition  entry  pressure  of  the  surrounding  media. 

This  analysis  leads  to  a  simple  expression  for  the  maximum  mound  thickness: 


T 


max  — 


Pe*  ~  Pe„nb 

kpg 


(1) 


where  Ped  is  the  water-NAPL  drainage  entry  pressure  of  the  aquifer  media,  Peimb  is  the 
water-NAPL  imbibition  entry  pressure  of  the  media,  Ap  is  the  density  difference  between 
the  DNAPL  and  water,  and  g  is  gravity.  This  is  similar  to  equations  given  by  Kueper  et 
al.  (1993)  and  McWhorter  and  Kueper  (1996)  with  two  differences  for  this  application. 
First,  they  use  the  entry  pressure  of  the  underlying  material  in  the  expression  to  give  the 
maximum  height  before  downward  DNAPL  penetration  into  the  aquitard.  Here  the  entry 
pressure  of  the  surrounding  aquifer  material  is  used  to  find  the  smaller  height  that  must  be 
reached  for  lateral  stability. 
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The  second  difference  in  these  expressions  is  the  capillary  pressure  assumed  for 
the  top  of  the  mound.  Kueper  et  al.  (1993)  and  McWhorter  and  Kueper  (1996)  both 
discuss  the  potential  water  imbibition  conditions  above  a  mound,  but  approximate  the 
capillary  pressure  as  close  to  zero  for  the  top  of  a  pool  height.  More  specifically, 
McWhorter  and  Kueper  (1996)  state  that  the  upper  surface  of  the  connected  DNAPL  lies 
where  the  capillary  pressure  is  zero.  That  is  probably  the  best  assumption  for  media  with 
a  wide  distribution  of  pore  sizes,  where  significant  connected  DNAPL  saturations  could 
remain  well  below  the  imbibition  entry  pressure,  and/or  where  the  media  might  not 
exhibit  a  distinct  imbibition  entry  pressure.  However,  in  media  with  a  uniform  narrow 
pore  size  distribution,  the  upper  boundary  of  high  DNAPL  saturation  can  only  correspond 
to  the  imbibition  entry  pressure.  That  is  the  value  that  should  be  used  in  the  expression  to 
estimate  the  maximum  static  thickness  on  a  horizontal  lens.  It  will  be  shown  later  that 
using  a  capillary  pressure  of  zero  for  the  top  of  the  mound,  instead  of  the  imbibition  entry 
pressure,  could  overestimate  the  thickness  by  2.5  times  in  certain  porous  media. 

The  mechanics  described  above  for  DNAPL  mounds  should  also  apply  to  static 
LNAPL  or  air  mounds  trapped  under  low  permeability  lenses  below  the  water  table. 
Figure  2.3  shows  the  hydrostatic  distribution  of  pressures  for  this  case.  Since  LNAPLs 
are  less  dense  than  water,  the  capillary  pressure  would  increase  with  elevation  to  a 
maximum  at  the  top.  Again,  the  advancing  edge  at  the  top  of  the  mound  would  be  at  the 
drainage  entry  pressure  of  the  media  while  the  bottom  receding  boundary  of  LNAPL 
saturations  correspond  to  pressures  on  the  water  imbibiton  curve. 
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Figure  2.3  Pressure  distribution  in  a  static  LNAPL  mound  trapped  below  the  water  table  beneath  a 
fine  textured  lens. 


Note  that  this  analysis  is  for  a  static  mound  on  a  hypothetical  flat  fine  soil  layer. 
It  does  not  explain  the  thickness,  and  thus  the  aquitard  penetrating  head,  that  could  build 
up  during  dynamic  flow  conditions  or  in  a  low  lying  depression.  Also,  it  only  predicts  a 
maximum  thickness,  assuming  the  NAPL  volume  and  infiltration  dynamics  were 
sufficient  to  reach  it.  Of  course,  a  flat  heterogeneity  is  not  necessary  for  this  analysis. 
The  mechanics  also  predict  the  additional  thickness  supported  above  and  in  addition  to 
DNAPL  pooled  in  a  depression.  Then,  this  analysis  would  predict  the  thickness 
supported  above  the  highest  point  of  a  depression  in  a  fine  lens. 

It  should  also  be  noted  that  this  analysis  is  for  a  material  with  a  non-zero  entry 
pressure,  i.e.,  one  that  exhibits  a  Brooks-Corey  Pc-S  relationship  (Corey,  1994). 
Macropores  may  exist  that  cannot  provide  the  necessary  capillary  resistance  for  a  static 
mound  of  significant  size.  A  van  Genuchten  Pc-S  relationship  (van  Genuchten,  1 980) 
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assumes  a  connected  NAPL  and  some  saturation  at  all  capillary  pressures.  However, 
even  if  a  van  Genuchten  porous  media  permits  some  flow  of  the  NAPL  below  the  entry 
pressure,  the  relative  permeability  of  the  NAPL  would  typically  be  so  low  that  the 
difference  in  mound  thickness  is  insignificant  or  only  manifests  itself  after  a  very  long 
time. 

This  analysis  has  implications  for  NAPL  remobilization  since  a  mound  at  a 
maximum  stable  thickness  would  clearly  be  at  the  point  of  incipient  motion.  First,  the 
entry  pressure  is  directly  proportional  to  the  interfacial  tension  between  the  two  fluids. 
Therefore,  a  lowering  of  the  interfacial  tension,  through  surfactant  flushing  for  example, 
could  cause  the  NAPL  to  again  become  mobile.  Second,  the  hydrostatics  show  that 
although  the  maximum  mound  thickness  is  independent  of  its  depth  below  the  water 
table,  changes  in  gradients  of  water  pressure  have  the  potential  to  remobilize  the  NAPL 
mound;  Indeed,  there  have  been  field  cases  in  which  DNAPL  that  had  been  static  for 
decades  mobilized  in  significant  volumes,  entering  wells  intended  only  for  pumping  and 
treating  the  aqueous  phase.  Groundwater  pumping  in  a  contaminated  zone  could  cause 
NAPL  migration  around  a  previously  perching  lens  and  deeper  into  the  aquifer.  The 
hydrostatics  also  show  a  possible  benefit  of  mass  removal.  Removing  NAPL  from  pools 
and  mounds  would  likely  reduce  their  potential  for  remobilization  and  redistribution. 


2.4  EXAMPLES  OF  THE  MAXIMUM  THICKNESS  OF  NAPL  MOUNDS 

The  above  analysis  illustrates  the  role  of  hysteresis  in  the  behavior  of  NAPLs  in 
porous  media,  and  how  that  behavior  enables  a  mound  of  significant  size  to  remain  static. 
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For  illustration,  table  2.1  contains  predicted  maximum  mound  thicknesses  for  various 
NAPLs  and  air  in  a  hypothetical  fine  and  medium  sand.  For  the  fine  sand,  an  air-water 
drainage  entry  pressure  of  70  cm  of  water  (6860  Pa)  and  an  air-water  imbibition  entry 
pressure  of  42  cm  of  water  (4116  Pa)  was  used.  For  the  medium  sand,  the  drainage  and 
imbibition  air- water  entry  pressures  were  assumed  to  be  20  cm  of  water  (1960  Pa)  and  12 
cm  of  water  (1 176  Pa),  respectively.  The  NAPL-water  entry  pressures  were  calculated  by 
scaling  the  air- water  entry  pressures  by  the  ratio  of  interfacial  tension  between  the  two 
fluids  (Lenhard  and  Parker,  1987).  Imbibition  entry  pressures  are  often  estimated  as  one- 
half  of  the  drainage  entry  pressure.  However,  in  the  experiments  presented  later  a 
Peimb/Ped  ratio  of  0.6  was  measured,  and  that  ratio  is  used  in  these  examples.  Note  that 
from  this  ratio  it  follows  that  if  Pc  =  0  is  assumed  for  the  upper  mound  boundary  instead 
of  Pc  =  Peimb,  equation  (1)  would  predict  a  maximum  mound  thickness  2.5  times  greater. 


Table  2.1.  Examples  of  the  maximum  thickness  of  laterally  stable  NAPL  mounds  in  medium  and  fine 
grained  sand  (air-water  entry  pressures  of  20  and  70  cm  of  water  (1960  and  6860  Pa),  respectively). 


Fluid 

Density  Interfacial 
Tension 

(g  cm-3)  (dynes  crn'l) 

Thickness  in 
Medium  Sand 

(cm) 

Thickness  in 
Fine  Sand 

(cm) 

Trichloroethylene 

1.47 

35 

8.3 

29.0 

T  etrachloroethy  lene 

1.63 

44 

7.8 

27.2 

1,1,1  -Trichloroethane 

1.35 

45 

14.3 

50.0 

1 ,2-Dichloroethane 

1.26 

30 

12.8 

44.9 

Air 

*0 

72 

8.0 

28.0 

Gasoline 

0.73 

50 

20.6 

72.0 

Benzene 

0.87 

35 

29.9 

104.7 

Jet  Fuel 

0.82 

50 

30.9 

108.0 

Creosote 

1.03 

20 

74.1 

259.3 
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2.5  LABORATORY  VERIFICATION 
Equipment 


Laboratory  experiments  were  performed  to  verify  equation  (1),  and  to  gain  insight 
into  the  configuration  the  mounds  would  exhibit.  Experiments  were  conducted  in  two 
thin  glass  flumes  (one  with  interior  dimensions  105  x  63  x  1  cm  and  the  second  with 
mtenor  dimensions  90  x  1 10  x  1  cm)  in  conjunction  with  a  light  attenuation  system. 
(Glass  et  al„  1989;  Tidwell  and  Glass,  1994).  Light  attenuation  measurements  take 
advantage  of  the  fact  that  the  intensity  of  light  transmitted  through  a  thin  slab  of  silica 
sand  is  strongly  affected  by  the  fluid  saturation.  The  thin  flumes  were  illuminated  from 
behind  and  images  captured  using  an  analog  ccd  camera  with  a  high  accuracy 
monochrome  frame  grabber.  This  method  allows  fluid  content  observations  to  be  based 
on  an  average  over  the  entire  thickness  of  the  experimental  cell,  reducing  the  influence  of 
occasionally  significant  edge  effects. 

Materials 

Two  sizes  (20/30  mesh  and  30/40  mesh)  of  Ottawa  sand  were  sieved  from  U.S. 
Silica  Sand  #14  (U.S.  Silica  Sand,  Ottawa,  IL)  and  used  in  the  experiments.  Sand  20/30 
was  used  in  experiments  1  and  3  and  sand  30/40  was  used  in  experiments  2  and  4.  The 
air-water  entry  pressures  for  the  sands  were  measured  in  the  same  experiment  flumes 
under  identical  packing  conditions  by  measuring  the  height  of  the  capillary  fringe  both  at 
equilibrium  drainage  and  imbibition  conditions.  These  are  given  in  table  2.2  along  with 
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the  hydraulic  conductivities  of  the  sands.  U.S.  Silica  Sand  #F95  (median  grain  size  150 
|xm)  was  used  for  the  fine  confining  layer  in  all  experiments. 

A  perfluorocarbon,  product  name  Performance  Fluid™  5070  (PF-5070, 3M 
Specialty  Chemicals,  St.  Paul,  MN),  was  used  for  the  DNAPL  in  experiments  1  and  2.  It 
was  selected  since  it  is  a  non-toxic  fluid  with  physical  properties  similar  to 
tetrachloroethylene  (described  in  Chapter  3).  Soltrol®  220  (Phillips  Petroleum  Co., 
Bartlesville,  OK),  a  hydrocarbon  commonly  used  in  LNAPL  experiments,  was  used  in 
experiments  3  and  4.  The  relevant  properties  of  these  fluids  are  given  in  table  2.3. 


Table  2.2.  Sand  Properties. 


20/ 30  Scmd 

30/40  Sand 

Ped  (cm  of  water) 

15.0 

21.0 

Peirob  (cm  of  water) 

8.9 

13.0 

Hydraulic  Conductivity  (cmr1) 

0.10 

0.06 

Table  2.3.  Fluid  Properties. 


PF-5070 

Soltrol  220 

Density  (g  cm*3) 

1.73 

0.79 

Viscosity  (centipoise) 

0.95 

3.65 

Interfacial  Tension  (dynes  cm-l)a 

a  .  *  .  . 

50.5 

35.2 

interfacial  tension  measured  with  fluids  dyed  and  at  an  equilibrium  value 
after  contact  with  water  in  the  experimental  flume.  The  water  surface  tension 
after  being  in  contact  with  the  sand  was  measured  to  be  approximately  68 
dynes  cm  -I. 


Experimental  Procedure 

Two  DNAPL  experiments  (1  and  2)  were  performed  where  DNAPL  was 
introduced  above  a  low  permeability  layer  from  a  well.  A  schematic  of  the  experiments 
is  shown  in  figure  2.4.  A  10  cm  fine  sand  layer  was  placed  at  the  bottom  of  the  flume. 
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Air  dry  sand  was  then  slowly  “rained”  in,  simultaneously  across  the  entire  flume  length, 
using  a  rectangular  funnel  and  coarse  screens  to  randomize  the  motion  of  the  sand  grains. 
This  method  achieved  a  very  uniform  packing  with  an  approximate  bulk  density  of  1.7  g 
cm 3.  C02  was  flushed  through  the  flume  to  displace  the  air  in  the  pore  space  in  order  to 
attain  a  high  water  saturation.  The  flume  was  then  slowly  filled  with  deionized  water 
dyed  with  FD&C  Red  #40  at  0.5  ml  per  100  ml  of  water  and  FD&C  Dyes  Yellow  #5  and 
Blue  #1  each  at  0.05  ml  per  100  ml  water.  The  DNAPL  used  in  these  experiments, 
PF-5070,  has  a  very  low  solvency  and  cannot  be  dyed.  Therefore,  this  dye  mixture  was 
necessary  to  minimize  light  transmission  through  the  water  phase  to  allow  observation  of 
the  DNAPL  mound. 

Wei! 


Figure  2.4.  Schematic  of  experiments  showing  how  NAPL  was  injected  through  the  well  and  the 
subsequent  redistribution  of  the  NAPL  mound. 
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A  well  ( 3  x  1  cm  rectangular  cross  section)  was  located  at  the  left  of  the  flume  for 
DNAPL  injection  and  to  monitor  the  capillary  pressure  during  the  experiment.  The 
PF-5070  was  injected  into  the  media  through  the  screened  well  by  maintaining  a  well 
thickness,  and  thus  a  capillary  pressure,  above  the  entry  pressure  of  the  media.  This 
method  allowed  a  mound  height  exceeding  the  predicted  maximum  thickness  to  be 
injected  while  minimizing  lateral  spreading  across  the  limited  length  of  the  experimental 
flume.  By  injecting  at  a  boundary  this  method  also  produces  a  line  of  symmetry, 
enabling  one  half  of  a  larger  symmetric  mound  to  be  observed  in  the  small  cell. 

Once  a  mound  of  sufficient  height  was  in  place,  injection  was  slowed  by 
decreasing  the  DNAPL  head  in  the  well,  but  still  leaving  it  above  the  predicted  final  well 
thickness  to  allow  it  to  equilibrate  on  its  own.  The  capillary  pressure  and  the  mound 
height  were  then  monitored  for  several  days,  and  in  one  case  for  several  weeks,  to  assure 
that  the  mound  had  assumed  a  stable  configuration. 

The  above  procedure  was  repeated  in  experiments  3  and  4  with  Soltrol  220,  dyed 
with  0.5  %  by  weight  Sudan  IV  organic  dye  (Aldrich  Chemical  Co.,  Milwaukee,  WI).  In 
these  experiments,  a  confining  layer  was  placed  across  the  top  of  the  media.  Similar  to 
the  procedure  in  the  DNAPL  experiment,  a  slug  of  the  LNAPL  was  injected  into  the 
media  along  the  well  and  under  the  confining  layer.  The  capillary  pressures  and  inverted 
mound  height  were  monitored  over  several  days  and  exhibited  behavior  similar  to  that  of 
the  DNAPL  experiments. 
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2.6  RESULTS 


The  measured  and  predicted  maximum  mound  heights  for  all  experiments  are 
given  in  table  2.4,  along  with  the  equilibrium  capillary  pressures  at  the  leading  edge  of 
the  mounds.  Figure  2.5  shows  this  same  information  on  a  graph.  The  values  predicted 
using  equation  (1)  correspond  well  to  the  measured  values  in  all  experiments.  The  error 
between  predicted  and  measured  values  is  well  within  the  measurement  error  of  the  input 
parameters.  Figures  2.6  and  2.7  are  photographs  of  the  NAPL  mounds  in  experiments  2 
and  4  respectively. 

The  NAPL  mounds  in  these  experiments  all  approached  a  static  equilibrium 
relatively  quickly.  Figure  2.8  shows  the  predicted  and  measured  DNAPL  mound  height 
with  time  along  with  the  approximate  capillary  pressure  at  the  mound  leading  edge  for 
experiment  1.  Experiment  1  shows  an  additional  DNAPL  injection  at  50  hours.  This  was 
done  to  assure  that  the  maximum  mound  thickness  had  been  exceeded.  The  DNAPL 
level  in  the  well  was  used  to  determine  the  capillary  pressure  at  the  bottom  leading  edge 
of  the  mound.  This  calculation  was  made  assuming  the  mound  was  at  equilibrium,  i.e.,  no 
horizontal  pressure  gradients.  These  measurements  confirmed  that  the  mound  did  stop  at 
the  oil-water  drainage  entry  pressure  of  the  media.  The  LNAPL  mounds  took  slightly 
longer  to  reach  a  static  equilibrium,  due  to  a  higher  viscosity  and  smaller  density 
difference,  but  exhibited  behavior  similar  to  that  of  the  DNAPL  experiments.  These  plots 
for  experiments  2,  3,  and  4  are  included  in  Appendix  B. 


11-32 


Table  2.4.  Maximum  NAPL  mound  thickness  in  laboratory  experiments 
compared  to  values  predicted  assuming  both  the  capillary  pressure  (Pc) 
at  the  top  of  the  mound  to  be  equal  to  zero  and  to  the  imbibition  entry 
pressure  (Peimb).  Also  shown  is  the  capillary  pressure  at  the  leading  edge 
of  the  mound  compared  to  the  predicted  oil-water  entry  pressure  of  the 
media. 


Predicted 
using  Pc  ~  0 

Predicted 
using  Pc  =  Peimb 

Measured 

Mound  Thickness 
(cm) 

Exp.  1 

15.2 

6.2 

8.0 

Exp.  2 

21.4 

8.1 

7.5 

Exp.  3 

37.0 

15.0 

16.5 

Exp.  4 

51.8 

20.3 

17.5 

Equil.  Cap.  Press. 

(cm  of  water) 

Exp.  1 

11.1 

11.4 

Exp.  2 

15.6 

14.3 

Exp.  3 

7.8 

8.4 

Exp.  4 

1 

10.8 

11.5 

Measured  Thickness  (cm) 


Figure  2.5.  Plot  comparing  measured  maximum  mound  thicknesses  with  predictions  based  both  on 
assuming  the  top  of  the  mound  occurs  at  a  capillary  pressure  equal  to  zero  and  equal  to  the  imbibition  entry 
pressure. 
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Figure  2.7.  Photograph  of  LNAPL  mound  (dark  area)  in  experiment  4  after  120  hours. 
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Pc  (cm  of  water)  Mound  Height  (cm) 


40 


Time  (hours) 


Figure  2.8.  Predicted  and  measured  DNAPL  mound  height  with  time  along  with  the  capillary  pressure  at 
the  mound  leading  edge  for  experiment  1. 
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Although  equilibrium  was  approached  rapidly  and  a  stable  mound  was  created  in 
these  experiments,  it  should  be  noted  that  this  was  not  the  case  in  experiments  where  the 
sand  contained  small  cracks  or  heterogeneities  of  coarse  sand.  It  these  cases,  the  mound 
still  rapidly  approached  the  predicted  height,  but  continued  to  slowly  decrease  as  the 
NAPL  moved  along  and  filled  the  crack.  If  this  pathway  of  low  entry  pressure  led  to  the 
opposite  well,  the  media  never  had  the  capillary  resistance  to  support  a  significant 
pressure,  and  the  mound  height  continued  to  diminish  as  it  filled  the  well.  This  would 
also  occur  in  homogeneous  media  when  the  advancing  mound  intersected  the  opposite 
well  due  to  having  too  much  volume  injected.  Once  this  pathway  of  NAPL  had  been 
established  across  the  media  boundary,  the  NAPL  always  seemed  to  be  able  to  flow  into 
the  opposite  well,  even  after  the  mound  reached  and  fell  below  the  entry  pressure  of  the 
media.  One  experiment  involving  air  was  attempted  and  also  showed  these  effects.  The 
air  did  mound  under  the  lens,  but  easily  found  or  even  formed  cracks  in  the  sand  that 
allowed  its  escape.  The  significance  of  these  effects  were  dependent  mainly  on  the  spill 
volumes  and  limited  size  of  the  experimental  cells.  However,  they  do  show  that  this 
theory  is  sensitive  to  minor  heterogeneities,  and  that  the  actual  static  height  reached  is 
dependent  on  the  volume  of  the  NAPL  spill  and  the  extent  of  spreading  through 
macroscale  pathways. 

These  effects  also  highlight  that,  unlike  in  these  controlled  homogeneous 
experiments,  NAPL  mounds  in  the  field  would  have  irregular  shapes.  The  mounds  did 
assume  similar  configurations  in  all  experiments.  The  maximum  thickness  of  the  mound 
occurred  only  near  the  point  of  injection,  with  the  mound  advancing  for  some  distance  in 
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a  narrow  finger.  This  shows  that  using  the  maximum  thickness  to  estimate  the  mound 
volume  might  have  limited  use  even  in  homogeneous  media. 

2.7  CONCLUSIONS 

Stable  mounds  of  NAPL  can  form  on  heterogeneities,  with  lateral  stability 
existing  even  on  flat  or  sloped  lenses.  Fine  textured  media  or  small  density  differences 
between  the  NAPL  and  water  can  lead  to  particularly  large  mounds.  This  behavior 
applies  both  to  DNAPL  mounds  perched  on  a  lens  or  the  bottom  of  an  aquifer  and  to 
LNAPLs  and  air  trapped  below  the  water  table  beneath  fine  textured  lenses. 

Modifying  previous  analyses  given  for  NAPL  mounds,  it  was  shown  that  the 
upper  boundary  of  the  mound  (the  boundary  of  high  NAPL  saturations)  in  uniform  media 
corresponds  more  closely  to  a  capillary  pressure  equal  to  the  imbibition  entry  pressure  of 
the  media  than  to  a  capillary  pressure  of  zero.  A  predictive  expression  for  the  maximum 
thickness  of  a  laterally  stable  NAPL  mound  was  given  and  verified  with  laboratory 
experiments.  This  immiscible  fluid  behavior  cannot  be  explained  without  considering  the 
hysteresis  in  capillary  pressure-saturation  relationships.  Models  not  incorporating 
hysteresis  would  likely  miss  this  phenomena  and  may  provide  incorrect  distributions  of 
NAPLs  in  the  subsurface. 
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Chapter  3 


A  Non-Toxic  DNAPL  Surrogate  for  Laboratory  and  Field  Experiments 

Abstract 

Laboratory  experiments  using  dense  chlorinated  solvents  present  concerns  with 
health,  disposal,  and  reactivity  with  laboratory  materials.  Regulatory  restrictions  prohibit 
field  experiments  with  chlorinated  solvents.  Perfluorocarbons  are  inert,  non-toxic, 
denser-than-water,  nonaqueous  phase  liquids,  most  having  viscosities  less  than  that  of 
water.  The  fluid  properties  of  this  family  of  chemicals  closely  match  those  of  chlorinated 
solvents.  The  properties  of  these  compounds  relevant  to  subsurface  transport  and  related 
research  were  determined  and  are  presented  here.  A  laboratory  spill  of  a  perfluorocarbon 
into  saturated,  heterogeneous,  porous  media  was  conducted  to  provide  qualitative 
observation  of  its  behavior  and  to  evaluate  its  visual  detection  during  an  experiment. 
Their  properties  show  that  perfluorocarbons  may  be  suitable  surrogate  fluids  for  many 
laboratory  and  field  investigations  into  the  behavior  of  dense  chlorinated  solvents  in 
groundwater. 

3.1  INTRODUCTION 

The  widespread  use  and  physical  properties  of  dense  chlorinated  solvents  have 
made  them  common  and  uniquely  problematic  groundwater  contaminants  (Pankow  et  al., 
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1996).  Much  research  is  currently  conducted  into  the  behavior  and  remediation  of  this 
group  of  denser-than-water,  nonaqueous  phase  liquids  (DNAPLs).  However,  chlorinated 
solvents  are  toxic,  volatile,  strong  solvents,  and  regulated  hazardous  materials. 

Therefore,  laboratory  experiments  with  them  present  concerns  for  human  health  and 
safety,  reactivity  with  laboratory  materials,  and  proper  disposal.  These  factors  often 
complicate  or  prohibit  needed  laboratory  experiments  (Kueper  and  Gerhard,  1995; 
Thomson  et  al.,  1992),  and  consequently,  researchers  have  commonly  desired  and  long 
sought  after  a  suitable  surrogate  fluid.  A  non-toxic,  non-reactive  DNAPL  would  benefit 
laboratory  research. 

Since  chlorinated  solvents  are  hazardous  materials,  regulations  prohibit  field 
experiments  with  them  (Kueper  and  Frind,  1991).  Except  for  a  number  of  notable 
experiments  in  the  Borden  aquifer  in  Canada  (e.g.,  Poulsen  and  Kueper,  1992;  Kueper  et 
al.,  1993)  most  research  has  been  limited  to  laboratory  experiments,  case  studies,  and 
numerical  simulations.  The  ability  to  conduct  field  experiments  with  a  non-toxic  DNAPL 
would  be  a  tremendous  asset  to  researchers. 

Perfluorocarbons,  a  non-toxic  family  of  compounds  with  densities  and  viscosities 
very  similar  to  those  of  dense  chlorinated  solvents,  are  identified  here  for  possible  use  as 
DNAPL  surrogates  in  laboratory  and  field  experiments.  Their  properties  most  relevant  to 
DNAPL  research  were  determined  and  are  presented. 
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3.2  PERFLUOROCARBONS 


The  perfluorocarbons  investigated  here  are  made  by  3M,  St.  Paul,  MN,  with  the 
product  names  of  Fluorinert™  liquids  and  Performance  Fluids™.  Fluorinert  and 
Performance  Fluids  are  chemically  identical  compounds,  except  that  the  Fluorinert 
liquids  are  subjected  to  higher  quality  assurance  for  use  in  different  applications  (personal 
communication,  Steve  Pignato,  3M  Specialty  Chemicals  Division,  St.  Paul,  MN,  1996). 

Fluid  Properties  of  DNAPLs.  DNAPL  flow  in  the  subsurface  is  governed  largely 
by  the  fluid  properties  density,  viscosity,  and  interfacial  tension.  The  high  densities  and 
low  viscosities  of  chlorinated  solvents  allow  relatively  rapid  and  extensive  downward 
movement  in  the  subsurface  (Pankow  et  al.,  1996).  DNAPLs  are  nonwetting  fluids  with 
respect  to  water  for  most  geologic  materials.  The  interfacial  tensions  between  chlorinated 
solvents  and  water  are  low  enough  to  permit  entry  into  small  fractures  and  pores,  but  at 
the  same  time,  high  enough  that  transport  is  strongly  affected  by  subsurface 
heterogeneities  (Pankow  et  al.,  1996).  These  properties  of  chlorinated  solvents  are 
closely  matched  by  various  perfluorocarbons.  Properties  of  some  common  DNAPLs  are 
given  in  table  3.1.  Table  3.2  lists  these  properties  for  various  perfluorocarbons. 


Table  3.1.  Properties  of  common  DNAPLs:  Trichloroethylene  (TCE),  Tetrachloroethylene  (PCE), 
and  1,1,1-Trichloroethane  (TCA).  (from  Cohen  and  Mercer,  1993) 


Density 

(g  cm'3) 

Absolute 

Viscosity 

(cP) 

“/*  - - - 

Kinematic 

Viscosity 

(cSt) 

Surface 
Tension 
(dynes  cm~l) 

Interfacial 
Tension 
(dynes  cm~l) 

Solubility 

In  Water 
(mgl-1) 

Vapor 
Pressure 
@~25  V 
(mm  Hg) 

TCE 

1.47 

0.57 

0.38 

29 

34 

1100 

57.9 

PCE 

1.63 

0.93 

0.57 

32 

44 

150 

17.8 

TCA 

1.35 

0.91 

0.67 

25 

45 

757 

123.0 
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Table  3.2.  Properties  of  various  3M  Fluorinert™  Liquids  (FC)  and  Performance  Fluids™ 


Product a 

Density 

(gem-3) 

Absolute 

Viscosity 

(cP) 

Kinematic 

Viscosity 

(cSt) 

Surface 
Tension 
(dynes  cm~l) 

Interfacial  b 
Tension 
(dynes  cm~l) 

Solubility 

In  Water 
(mg  l-l) 

Vapor  c 
Pressure 
(mm  Hg) 

FC-75 

1.80 

1.40 

0.77 

15 

<10 

31 

FC-87 

1.63 

0.65 

0.40 

9 

<10 

610 

FC-72 

1.68 

0.67 

0.40 

12 

57 

<10 

232 

FC-84 

1.73 

0.95 

0.55 

13 

58 

<10 

79 

FC-40 

1.87 

4.11 

2.2 

16 

52 

<10 

3 

PF-5060d 

1.68 

0.67 

0.40 

12 

<10 

232 

PF-5070' 

1.73 

0.95 

0.55 

13 

56 

<10 

79 

PF-5080f 

S  T'liinri I  \ 

1.76 

nmdc  oro 

1.40 

t  ntarl  Im  r  17/” 

0.77 

■»  n — - — 

15 

n..:  i _ r»r? 

<10 

45 

- - - - v  V  J  A  Wj  M.  VilVUIIIMlVV  1  IU1UJ  VJ  J  X  A  m 

b  measured  at  18°  C  using  a  ring  tensiometer  (CSC  Scientific  Company,  Fairfax,  VA) 
c  at  25°  C 
d  same  as  FC-72. 
e  same  as  FC-84. 

f  same  as  FC-75,  although  values  differed  slightly  between  sources. 


Chemical  Properties.  Fluorinert  liquids  and  Performance  Fluids  are  clear, 
colorless,  and  odorless.  They  consist  primarily  of  fluorinated  chains  of  six  to  eight 
carbons.  Fluorination  is  complete,  so  they  contain  no  hydrogen  or  chlorine,  and  therefore 
do  not  contribute  to  stratospheric  ozone  layer  damage  (3M,  1994).  The  very  strong 
carbon- fluorine  bonds  result  in  these  compounds  being  chemically  inert,  and  extremely 
stable  (Ju  et  al.,  1991).  They  are  non-flammable,  and  will  exert  essentially  no  solvent 
action  on  non-fluorinated  compounds  (3M,  1994).  They  have  very  low  water  solubilities, 
and  are  immiscible  with  most  organic  compounds  as  well  (Zhu,  1993).  One  particularly 
interesting  characteristic  is  that  these  perfluorocarbons  have  high  oxygen  solubilities,  on 
the  order  of  50  ml  of  gas  per  100  ml  of  liquid  at  1  atm  and  25°C  (Ju  et  al.,  1991). 
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Health  and  Safety.  The  Material  Safety  and  Data  Sheet  (MSDS)  for  Performance 
Fluid  5070  (PF-5070)  states  that  this  fluid  poses  no  significant  health  hazard  unless 
exposed  to  extreme  heat.  At  high  temperatures,  some  toxic  decomposition  products  can 
occur.  As  a  “good  industrial  hygiene  practice”,  the  MSDS  recommends  proper 
ventilation,  but  no  irritation  is  expected  from  moderate  inhalation.  PF-5070  is  listed  as 
non-irritating  for  dermal  and  eye  exposure,  and  “practically  non-toxic”  for  oral  exposure. 
PF-5070  poses  insignificant  toxicity  to  aquatic  organisms,  including  waste  treatment 
microorganisms. 

The  medical  uses  of  similar  perfluorocarbons  give  further  testimony  to  their  low- 
toxicity.  Their  high  oxygen  solubility,  combined  with  being  inert  and  non-toxic,  give  rise 
to  applications  such  as  liquid  ventilation  for  human  lung  disease,  artificial  blood,  and  as 
an  oxygen  transfer  mechanism  in  bioreactors  (Goldfarb,  1992;  Leach  et  al.,  1996;  Spence, 
1995;  Ju  et  al.,  1991).  Perfluorocarbons  also  have  applications  in  electronics 
manufacturing  and  as  alternatives  to  chlorofluorocarbons  in  refrigeration. 

Permission  was  obtained  from  the  local  wastewater  treatment  plant  for  disposal  of 
small  quantities  in  the  municipal  sewer  system  (personal  communication,  David  Meyer, 
City  of  Fort  Collins  Wastewater  Treatment  Plant,  1996).  Furthermore,  preliminary 
permission  was  obtained  to  conduct  controlled  spills  of  PF-5070  in  field  experiments  in 
Colorado  (personal  communication,  Peter  Laux,  Colorado  Department  of  Public  Health 
and  Environment,  1996).  Of  course,  permission  would  have  to  be  obtained  on  an 
individual  basis  and  from  the  appropriate  local  agencies  for  experiments  and  disposal 
elsewhere. 
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Detection.  As  a  result  of  being  highly  inert,  perfluorocarbons  cannot  be  dyed. 

The  common  organic  dyes  Sudan  III  and  Sudan  IV  (Aldrich  Chemical,  Milwaukee,  WI) 
as  well  as  Oil  Blue  N  (Aldrich  Chemical)  and  Uvitex-OB  (Ciba-Geigy),  were  all  tested. 

It  was  determined  that  no  dye  will  be  dissolved  in  these  perfluorocarbons  except  for  a 
highly  toxic  dye  illegal  in  the  U.S.  (personal  communication,  Steve  Pignato,  3M,  St.  Paul, 
MN,  1995).  Visual  detection  in  the  laboratory  was  accomplished  by  using  dyed  water 
along  with  a  clear  Performance  Fluid. 

The  relative  dielectric  strength  of  the  Performance  Fluids  (around  1.8)  makes 
them  suitable  for  subsurface  monitoring  with  ground  penetrating  radar  and  time  domain 
reflectometry  as  discussed  by  Brewster  et  al.  (1995)  and  Kueper  et  al.  (1993).  Also,  the 
absence  of  perfluorocarbons  in  the  natural  environment  have  made  them  excellent 
chemical  tracers  in  applications  such  as  monitoring  building  airflow  and  in  studies  of 

borehole  cross-contamination  in  groundwater  wells  (Cheong  and  Riffat,  1995;  Russell,  et 
al.,  1989). 

Gamma  Attenuation.  Gamma  energy  attenuation  has  been  used  as  a  laboratory 
measurement  technique  in  a  number  of  studies  of  DNAPL  flow  and  distribution  in  soil 
(e.g.,  Ferrand  et  al.,  1986;  Illangasekare  et  al.,  1995).  By  knowing  how  much  a  fluid 
attenuates  gamma  energy,  the  fluid  saturations  in  a  soil  column  can  be  determined.  The 
cesium  (  Cs)  and  americium  (24lAm)  gamma  attenuation  coefficients  for  a  few 
perfluorocarbons  were  measured  and  are  given  in  table  3.3  along  with  the  attenuation 
coefficients  for  water,  trichloroethylene,  and  1,1,1-trichloroethane. 
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Table  3.3.  Gamma  Attenuation  Coefficients,  (cm2  g'1) 


137Cs 

241  Am 

FC-84  * 

0.073 

0.174 

FC-72 

0.076 

0.172 

1,1,1  Trichloroethane  (TCA) b 

0.078 

0.336 

Trichloroethylene  (TCE) c 

0.070 

0.367 

Water 

0.092 

0.205 

a  FC  =  Fluorinert  liquid,  FC-84  same  product  as  PF-5070 
bfrom  Illangasekare  et  al.,  1995 
c  from  Ferrand  et  al.,  1986 


Water  Solubility.  Perfluorocarbons  have  very  low  solubilities  in  water  (<10 
mg/1).  This  is  a  limitation  for  some  research  applications,  such  as  in  monitoring 
dissolved  phase  contaminants  or  where  dissolution  of  a  NAPL  is  under  investigation. 
However,  this  is  an  advantage  in  other  applications.  For  example,  the  interfacial  tension 
between  water  and  a  NAPL  often  changes  with  time  as  the  dissolved  concentrations 
equilibrate,  often  complicating  analysis  (e.g.,  Schroth  et  al.,  1995).  The  interfacial 
tension  between  water  and  various  perfluorocarbons  was  monitored  over  a  period  of  7 
days  using  a  ring  tensiometer  (DuNouy  Tensiometers,  CSC  Scientific  Company,  Inc., 
Fairfax,  VA).  The  interfacial  tension  did  not  change  significantly  with  time,  decreasing 
less  than  4%. 

Cost.  The  Performance  Fluids  currently  cost  about  $50  per  liter.  This  compares 
with  about  $8  per  liter  for  trichloroethylene  (TCE)  and  about  $6  per  liter  for 
tetrachloroethylene  (PCE).  However,  the  disposal  costs  for  chlorinated  solvents  can 
range  from  $3  per  liter  to  $18  per  liter,  with  little  or  no  disposal  costs  for 
perfluorocarbons  if  sewer  disposal  is  permissible.  Also,  since  the  perfluorocarbons  are 
highly  inert  and  have  low  solubilities  of  water,  contamination  of  the  fluid  is  not  a  strong 
concern,  making  recovery  and  reuse  possible  in  the  laboratory. 
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3.3  EXPERIMENTAL  OBSERVATIONS 


A  laboratory  spill  of  Performance  Fluid  5070  (PF-5070)  into  a  saturated, 
heterogeneous,  sand  column  was  conducted.  This  was  done  primarily  to  evaluate  how 
well  the  clear  DNAPL  could  be  visually  detected  by  dying  the  water  phase  instead  of  the 
NAPL  as  is  typically  done,  and  also  to  gain  general  observations  and  experience  with  its 
use  in  laboratory  experiments.  PF-5070  was  selected  for  experimentation  since  its 
density  and  viscosity  most  closely  match  that  of  PCE. 

The  experiment  was  conducted  in  a  thin  glass  flume  (1  cm  x  63  cm  x  107  cm)  in 
conjunction  with  a  light  attenuation  system  similar  to  that  described  by  Glass  et  al.  (1989) 
and  Tidwell  and  Glass  (1994).  The  column  was  packed  with  a  medium  sand  (median 
grain  diameter,  D50  =  0.50  mm)  containing  two  very  coarse  sand  lenses  (D50  =  2.00  mm) 

and  a  fine  sand  lens  (D50  =  0.15  mm).  A  schematic  of  the  sand  packing  is  shown  in  figure 
3.1. 


Figure  3.1.  Schematic  of  sand  packing  for  experiment  with  Performance  Fluid  5070. 
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Light  attenuation  measurements  take  advantage  of  the  fact  that  the  intensity  of 
light  transmitted  through  a  thin  column  of  silica  sand  is  strongly  affected  by  the  fluid 
saturation.  Since  both  the  clear  DNAPL  and  water  easily  transmit  light  through  the 
column,  the  water  was  dyed  with  red  food  color  (FD&C  Red  #40)  to  a  concentration  of 
0.5  ml  per  100  ml  water  plus  green  food  color  (FD&C  Dyes  Yellow  #5  and  Blue  #1)  at 
0.05  ml  per  100  ml  water.  The  purpose  of  this  mixture  was  to  minimize  light 
transmission  through  the  water  phase,  and  therefore,  allow  visual  DNAPL  monitoring.  In 
addition  to  color  photographs,  images  were  captured  using  an  analog  ccd  camera  (Cohu 
4910  Series,  San  Diego,  CA)  in  conjunction  with  a  high  accuracy  monochrome  frame 
grabber  and  image  analysis  software  (Data  Translation,  Inc.,  Model  DT3155  frame 
grabber  and  Global  Lab  Image™  software,  Marlboro,  MA). 

Figure  3.2(a)  shows  the  saturated  column  before  the  DNAPL  spill.  Figures  3.2(b) 
and  3.3  were  taken  during  the  spill.  The  bright  areas  in  figure  3.2(b)  indicate  the  presence 
of  the  clear  DNAPL.  The  red  areas  are  the  dyed  water.  As  expected,  the  perfluorocarbon 
fluid  required  some  positive  capillary  pressure  to  penetrate  the  saturated  zone  and  then 
moved  easily  downward.  In  the  saturated  zone,  the  coarse  sand  was  a  preferential  flow 
path  and  trap  to  the  fluid,  while  the  fine  lens  acted  as  a  barrier.  This  is  behavior  typical  of 
DNAPLs  (Illangasekare  et  al.,  1995). 

In  the  coarse  material,  the  presence  of  the  DNAPL  was  very  clear.  However,  the 
medium  sand  had  a  higher  water  content,  causing  a  less  distinct  difference  between  the 
light  transmitted  by  the  two  fluids,  and  therefore,  making  DNAPL  detection  more 
difficult.  Furthermore,  some  edge  effects  were  observed  in  the  medium  sand,  i.e., 
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Figure  3.2.  Experimental  column  before  (a)  and  after  (b)  the  spill  of  PF-5070. 
The  water  is  dyed  red  while  the  bright  areas  indicate  high  DNAPL  saturations. 


DNAPL  preferentially  moving  along  the  front  and  rear  glass.  One  advantage  of  using  a 
thin  column  and  light  attenuation  measurements  is  that  it  allows  observation  of  the  fluid 
saturations  averaged  over  the  column  thickness.  This  does  allow  more  accurate 
determination  of  the  significance  of  any  edge  effects.  Although  more  difficult  in  this 
medium  sand,  DNAPL  detection  was  still  possible.  Figure  3.3  is  a  close  up  taken  with 

the  image  analysis  equipment  of  the  DNAPL  in  the  medium  sand  as  it  mounded  on  and 
flowed  around  the  fine  lens. 


Figure  3.3.  A  close-up  view  of  the  laboratory  spill, 

it  mounded  on  and  flowed  around  a  fine  lens. 


showing  the  DNAPL  presence  in  the 


medium  sand  as 


3.4  CONCLUSIONS 


A  group  of  non-toxic,  inert,  dense  nonaqueous  phase  liquids  has  been  identified 
and  investigated  as  sutrogate  fluids  for  chlorinated  so, vents  in  laboratory  and  field 

experiments.  They  have  fluid  propetties  very  similar  to  those  of  flte  dense  chlorinated 
solvents  bu,  pose  no  heahh  hazard  and  are  no,  reguIa,ed  hazardous  materials.  They 
no,  reactive  with  laboratoty  materials.  Therefote,  these  perfluorocarbon; 
asset  to  laboratoty  and  field  research  on  dense  chlorinated 


are 

s  could  be  a  great 


solvents  in  groundwater. 


Note:  An  additional  fluorocarbon  which  only  recently  became 


mentioned  in  Appendix  C  along  with 


commercially  available  is 


some  ideas  for  further  research  on  these  fluids.) 
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Chapter  4 


GEOMETRY  OF  LNAPL  MOUNDS  IN  WATER-NAPL-AIR  SYSTEMS 
Abstract 

Light  nonaqueous  phase  liquids  (LNAPLs),  such  as  petroleum  hydrocarbons,  are 
released  to  the  subsurface  environment  through  spills  or  leaks  from  buried  storage  and 
transmission  facilities.  When  spillled  at  or  near  the  land  surface,  an  LNAPL  will  migrate 
vertically  downward  with  some  lateral  spreading,  until  it  reaches  a  water-saturated  zone. 
If  sufficient  volume  of  LNAPL  has  been  spilled,  it  will  form  a  lens  and  spread.  This  paper 
presents  a  conceptual  model  of  an  LNAPL  lens  and  summarizes  laboratory  experiments 
that  illustrate  the  equilibrium  distribution  of  an  LNAPL  lens  at  the  capillary  fringe  in  a 
homogeneous  sand. 

4.1  INTRODUCTION 

The  contamination  of  soils  and  aquifers  by  LNAPLs  is  a  widespread  problem. 
Gasoline  leaked  from  underground  storage  tanks  is  an  example  of  an  LNAPL  that  poses 
potential  environmental  risks  in  both  the  vadose  zone  and  underlying  saturated  zone. 
Predicting  the  movement  and  distribution  of  LNAPLs  in  the  subsurface  environment  is 
critical  to  the  design  of  effective  mitigation  measures  for  contaminated  soils. 

Capillary  pressure-saturation  relationships  in  a  porous  medium  determine  the 
distribution  of  immiscible  fluids  under  static  conditions  and  can  largely  influence  the 
movement  of  the  fluids  when  the  system  is  not  at  equilibrium.  Capillary  pressure  is  a 
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function  of  the  physical  and  chemical  properties  of  the  fluids  and  the  porous  media,  and  is 
defined  as: 

Pc  =  Pnw  -  Pw=  2a  cosa/r  (4_1) 

where  Pc  is  the  capillary  pressure,  P„w  is  the  pressure  in  the  nonwetting  phase,  Pw  is  the 
pressure  in  the  wetting  phase,  a  is  the  interfacial  tension  between  the  fluids  in  contact;  a  is 
the  contact  angle  of  the  fluid-fluid  interface  with  the  porous  media,  and  r  is  the  effective 
radius  of  curvature  of  the  interface  (Corey,  1994).  In  a  three-phase  water-NAPL-air 
system,  it  is  usually  assumed  that  the  water  phase  is  the  wetting  phase,  the  NAPL  phase  is 
intermediate,  and  the  air  phase  is  the  nonwetting  phase  (Lenhard  and  Parker,  1987).  In 
the  presence  of  three  fluids,  two  relevant  phase  interfaces  exist:  a  water-NAPL  interface 
and  aNAPL-air  interface  (van  Geel  et  al.,  1994). 

Capillary  pressure-saturation  (Pc-S)  relationships  are  hysteretic,  i.e.  fluid 
saturations  or  contents  are  not  unique  functions  of  the  capillary  pressure.  Several  studies 
have  shown  the  importance  of  hysteresis  in  multifluid  flow  in  porous  media  (Lenhard, 
1992;  van  Geel  and  Sykes,  1994;  Essaid  et  al.,  1993;  and  Buesing  and  Pantazidou,  1996). 

Pantazidou  and  Sitar  (1993)  conducted  two-dimensional  LNAPL  infiltration 
experiments  in  sands  and  derived  equations  to  predict  the  maximum  vertical  dimension  of 
LNAPL  lenses  in  the  vadose  zone.  In  their  model,  they  simplified  the  actual  pore 
geometry  of  the  sand  by  assuming  that  NAPL  advance  and  drainage  was  dictated  by  the 
sizes  of  pore  necks  and  pore  bodies,  each  of  which  could  be  identified  with  a  characteristic 
diameter.  They  considered  a  lens  with  curved  upper  and  lower  boundaries.  To  compute 
the  vertical  dimensions,  Pantazidou  and  Sitar  (1993)  wrote  equations  for  the  balance  of 
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forces  acting  at  three  points  on  the  lens  boundaries.  The  midpoint  on  the  upper  boundary 
of  the  lens  was  assumed  to  be  an  air-NAPL  interface.  Because  NAPL  is  the  wetting  fluid 
with  respect  to  air  and  drains  from  the  pores,  this  interface  was  located  in  the  pore  necks. 


Governing  Forces  at: 
Point  1)  NAPL- Air 
Point  2)  NAPL-Water 
Point  3)  NAPL- Air 


1 


Figure  4-1.  Schematic  of  an  LNAPL  lens  at  equilibrium,  as  presented  by  Pantazidou  and  Sitar 
(1993). 

The  midpoint  on  the  lens  lower  boundary  was  assumed  to  be  governed  by  a  NAPL-water 
interface.  It  was  assumed  water  is  more  wetting  with  respect  to  the  NAPL  and  drains 
from  the  pores,  thus  this  interface  should  also  be  located  in  the  pore  necks.  The  difference 
in  NAPL  pressures  between  points  1  and  2  at  equilibrium  is  given  by  the  static  NAPL 
pressure  distribution.  Solving  for  T,  Pantazidou  and  Sitar  (1993)  derived  the  equation: 

T  =  (l/p0g)[4(Oow  +  oa0)/d„  -  pwghw]  (4-2) 
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where  T  is  the  maximum  vertical  thickness  of  the  NAPL  lens;  and  p0  is  the  NAPL  density, 
pw  is  the  water  density;  aow  and  aao  are  the  interfacial  oil-water  tension  and  oil  surface 
tension,  respectively,  and  hw  is  the  distance  between  the  water  table  and  the  bottom  of  the 
lens.  They  assumed  the  top  boundary  was  oil-air  drainage,  the  bottom  boundary  was 
water-oil  drainage  and  the  advancing  edge  of  the  lens  was  at  the  oil-air  imbibition  entry 


pressure. 


Schroth  et  al.  (1995)  extend  the  theory  of  Pantazidou  and  Sitar  (1993)  by  making 
the  assumption  that  the  top  of  the  lens  was  at  the  top  of  the  water-air  capillary  fringe.  The 
pore  neck  diameter  d„  was  calculated  using  the  Laplace  capillary  rise  equation: 


d„  =  (4aaw)/(Pwghcap) 


(4-3) 


where  h^p  is  the  capillary  fringe  height.  The  height  of  point  2  above  the  water  table,  h„ 
(Figure  4-1)  was  expressed  as: 

hw  =  heap  -  T  (4-4) 

Eq.  4-4  was  inserted  into  Eq.  4-2.  The  total  lens  thickness  then  becomes: 

T  =  l/p0g[4(aa0  +  cow)/dn  -  Pwgh^p  +  pwgT]  (4_5 


(4-5) 


Solving  Eq.  4-3  for  hcP  and  substituting  it  into  Eq.  4-5  and  solving  for  T  gave  a  folly 
explicit  equation  for  computing  T: 

T  —  4[(oaw  -  oa0  -  <j0w)]/[dng(pw  -  p0)]  (4_6 

Schroth  et  al.  (1995)  also  found,  in  their  experiments,  that  the  upper  boundary  was  flat. 
The  results  presented  by  Schroth  et  al.  (1995)  did  not  match  predictions  well.  They 
attributed  the  error  to  the  time  dependence  of  the  interfacial  tension. 
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In  contrast  to  Pantazidou  and  Sitar  (1993),  Abdul  (1988)  predicted  that  NAPLs 
form  lenses  in  the  capillary  fringe  with  horizontal  upper  and  curved  lower  boundaries.  He 
speculated  that  NAPL  would  fully  penetrate  the  capillary  fringe  and  at  equilibrium  form 
“pancake-shaped”  lenses. 

In  this  paper,  we  incorporate  three  phase  air-NAPL-water  relationships  in  a  new 
analysis  and  develop  a  conceptual  model  of  the  LNAPL  lens  during  infiltration  and  at 
equilibrium.  Two-dimensional  flume  experiments  were  used  to  verify  the  conceptual 
model. 


4.2  MATERIALS  AND  METHODS 

The  experimental  flume  (Figure  4-2)  was  constructed  from  two  sheets  of  122  cm  x  91cm  x 
1.27  cm  tempered  glass  placed  inside  of  a  steel  frame.  The  glass  was  separated  by  a  1-cm 
thick  aluminum  spacer,  which  was  sealed  to  the  glass  with  Silicone®  high  vacuum  grease. 
Two  vertical  bars,  the  height  of  the  column,  and  a  horizontal  bar  connecting  the  vertical 
bars,  were  used  to  minimize  deflection  of  the  glass  during  the  experiments.  Fully  screened 
vertical  wells  were  installed  at  each  end  of  the  flume.  A  light  box  was  set  up  behind  the 
flume  and  used  to  distinguish  the  saturated  and  unsaturated  zones.  An  analog  charge 
couple  device  (CCD)  camera  placed  in  front  of  the  flume  recorded  the  movement  of 
LNAPL  over  time.  A  framegrabber  converted  the  LNAPL  images  to  digital  format. 
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Figure  4-2.  Flume  Schematic  (not  to  scale). 

The  LNAPL  used  in  the  experiments  was  Soltrol  220®  (Phillips  Petroleum  Co., 
Bartlesville,  OK).  Soltrol  220®  is  a  mixture  of  branched  alkanes  having  low  solubility  in 
water  (Lenhard,  1992).  Its  components  can  be  found  in  a  variety  of  petrochemical 
products,  e.g.  diesel  fuel  (Schroth  et  al.,  1995).  To  distinguish  between  water  and  the 
LNAPL  during  the  experiments,  the  LNAPL  was  dyed  by  adding  0.01%  by  weight  of 
Sudan  IV®  (Aldrich  Chemical  Company,  Milwaukee,  WI).  Sudan  IV®  is  a  fat  stain  of 
bright  red  color  and  is  virtually  insoluble  in  water.  The  density  of  the  dyed  Soltrol  220® 
was  measured  in  the  laboratory  as  0.79  gm/cm3.  The  aqueous  phase  used  was  deionized 
water.  The  density  of  the  water  was  assumed  to  be  1  gm/cm3. 
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Two  sets  of  fluids  pairs  were  used  in  the  experiments.  Interfacial  tensions  are 
given  in  Table  4-1.  The  values  shown  were  measured  with  fluids  drained  from  the  flume 
after  the  experiments.  The  LNAPL-water  interfacial  tension  before  contact  with  the  sand 
was  44.6  dynes/cm. 

TABLE  4-1.  Values  of  Parameters. 


Fluid  Set  1  (Experiments  1-3) 

Parameter 

Air-water 

LNAPL- 

water 

Air-LNAPL 

LNAPL-water 

Units 

P«d 

20580 

10976 

8045 

5457 

P  1 

13720 

7317 

5363 

3638  ' 

Pci2 

14700 

7840 

5746 

3898 

a3 

66 

35.2 

25.8 

17.5 

| 

hed 

21 

11.2 

8.2 

5.6 

cm  of 

wetting  fluid 

hei* 

14 

7.5 

5.5 

3.7 

cm  of 

wetting  fluid 

h*i2 

15 

8.0 

5.9 

4 

cm  of 

wetting  fluid 

1.  Using  h«i  a'w  =  14  cm  estimated  in  flume. 

2.  Using  h«i  *'w  =  15  cm. 

3.  ASTM  Method 


A  spreading  coefficient  Csp,  is  defined  as: 

Csp  =  Oaw  -  (aso  +  CTow)  (4_7) 

where  a,w  ,a,0  and  ctow  are  the  surface  tensions  of  water  and  oil  and  the  oil  water 
interfacial  tension.  If  the  spreading  coefficient  is  positive,  the  oil  phase  will  spread  as  film 
flow  between  the  air  and  water  phases  in  the  unsaturated  zone.  The  spreading  coefficients 
are  10.3  and  28  for  fluid  sets  1  and  2,  respectively. 

Experiments  were  conducted  using  F-14  sand  (U.S.  Silica  Sand,  Ottawa,  Illinois). 
The  particle  size  distribution  was  measured  in  the  laboratory  (Figure  4-3)  by  sieving 
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(ASTM  Method  4464).  In  order  to  obtain  a  more  homogeneous  sand,  the  F-14  sand  was 
sieved  with  a  U.S.  Standard  #40  sieve  (0.425  mm).  Particles  less  than  0.425  mm  and 
greater  than  0.589  mm  were  discarded. 

A  2-cm  thick  layer  of  #8-12  (2.362-1.397  mm)  mesh  silica  sand  (Colorado  Silica 
Sand  Inc.,  Colorado  Springs,  Colorado)  was  placed  at  the  bottom  of  the  flume  to  facilitate 
drainage.  A  “rainer”  was  used  to  pack  the  flume  (Rad  and  Turney,  1987).  The  “rainer” 
consisted  of  an  acrylic  hopper  that  extended  across  the  full  width  of  the  flume  and  packing 
was  conducted  by  raining  continuously  to  prevent  local  gradations  in  particle  size.  The 
bulk  density  of  the  soil  was  approximately  1.7  gm/cm3. 

The  flume  was  initially  flushed  with  CO2  gas  through  the  bottom  ports  on  both 
sides  (Figure  4-2)  to  replace  the  air  in  the  pore  space,  and  then  slowly  saturated  with 
deionized,  deaired  water  from  the  bottom.  After  saturating  the  sand  pack,  the  water 
reservoir  was  lowered  to  a  level  of  8.5  cm  above  the  bottom  of  the  flume.  The  water  was 
allowed  to  drain  for  at  least  1 5  hours,  resulting  in  a  capillary  fringe  height  (h«j*'w  ) 
measured  visually  in  the  flume  of  21  cm  +  0.5  cm.  The  capillary  fringe  was  determined 
visually  as  the  zone  above  the  water  table  where  the  effective  saturation  is  very  close  to 
1 .0.  We  use  the  terms  capillary  fringe  and  entry  pressure  interchangeably.  Hence,  the 
entry  pressure  head  on  the  drainage  curve  for  the  sieved  sand  was  21  cm  ±  0.5  cm  (Pe<i1'w 
=  20580  dynes/cm2)  .  To  obtain  the  entry  pressure  head  on  the  imbibition  curve,  the  water 
level  was  raised  at  a  steady  rate  using  a  peristaltic  pump  to  a  level  of  40  cm  above  the 
bottom  of  the  flume.  The  water  was  allowed  to  imbibe  for  17  hours,  creating  a  fringe 
height  of  14  cm  +  0.5  cm. 


11-62 


Figure  4-3.  Particle  size  distribution  for  F-14  sand. 


4. 3  CONCEPTUAL  MODEL  OF  AN  LNAPL  SPILL  AND  LENS 

Figure  4  shows  the  ccd  camera  pictures  of  a  spill  with  fluid  set  2.  These  photos 
will  be  used  to  illustrate  the  conceptual  model.  When  an  LNAPL  is  spilled  near  the  soil 
surface,  it  moves  through  the  unsaturated  zone  as  a  result  of  gravity  and  pressure 
gradients,  (figure  4-4a).  As  the  LNAPL  reaches  areas  of  higher  water  saturation  near  the 
capillary  fringe,  it  spreads  horizontally  until  the  capillary  forces  equilibrate.  The  capillary 
pressure  at  the  advancing  edge  of  the  lens  is  the  water-oil  entry  pressure  on  the  drainage 
curve. 

An  LNAPL  with  a  positive  spreading  coefficient  spreads  horizontally  between  the 
water  and  air  phases  as  a  film  (Figure  4-4b).  The  LNAPL  continues  migrating  downward. 
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on  an  oil-water  drainage  curve  (Figure  4-4b,  point  A)  until  it  establishes  itself  on  the  top 
of  a  compressed  capillary  fringe  (Figure  4-4c,  point  B).  Over  time  the  LNAPL 


“rebounds”  upward  and  shifts  to  an  oil-water  imbibition  curve  at  the  bottom  of  the  lens 
(Figure  4c,  point  C),  as  it  continues  to  move  laterally  on  top  of  the  capillary  fringe,  along 
the  oil-water  drainage  curve  (Figure  4-4c,  point  D).  A  schematic  of  a  lens  at  equilibrium 
and  the  equilibrium  pressure  distribution  are  shown  in  figure  4-5. 


LNAPL 


Figure  4-4 a.  Injection  of  LNAPL. 


Figure  4-4b.  LNAPL  lcn  at  an  early  time. 
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Figure  4-4c.  LNAPL  lens  at  equilibrium 


Figure  4-5.  Pressure  distribution  and  notation  assumed  for  LNAPL  lens. 


4.4  LABORATORY  FLUME  EXPERIMENTS 


Two-fluid  phase  experiments  were  conducted  initially  to  verify  the  assumption  that 
the  advancing  edge  of  the  LNAPL  lens  was  at  a  capillary  pressure  equal  to  the  drainage 
entry  pressure  in  an  LNAPL-water  system.  The  flume  was  filled  with  the  sieved  sand  and 
saturated  with  deionized  water.  The  top  of  the  coarse  sand  was  at  a  height  of  60.5  cm 
above  the  bottom  of  the  flume.  A  layer  of  F-95  (U.S.  Silica  Sand,  Ottawa,  Illinois)  sand 
was  wet-packed  on  top  of  the  coarse  sand.  The  fine  sand  created  a  boundary  along  which 
the  LNAPL  would  advance  and  eventually  reach  equilibrium.  The  water  level  was  kept  at 
the  top  of  the  fine  layer  as  it  was  packed  in  the  flume.  The  wet  packing  was  used  to 
eliminate  any  entrapped  air  at  the  interface. 

Using  the  scaling  relationship  as  defined  by  Parker  et  al.  (1987),  the  entry  pressure 
head  for  the  LNAPL-water  drainage  curve,  hed°'w  was  estimated  as  1 1.2  cm  water  by: 

hed°'w  =  heda‘w  (3aw/0Ow)  (4-9) 

where  paw=l  and  P0w=OaW/cw  The  entry  pressure  head  for  the  LNAPL-water  imbibition 
curve,  hej°"w ,  was  estimated  as  3.7  cm  water  by: 

hei°"W  =  hei*  W  (Paw/pow)  (4-11) 

The  maximum  thickness  predicted  in  cm  of  oil  was  obtained  by: 

Tmax  =  hed°‘w  -  hej°"w  (in  cm  of  water)(pw/Ap)  (4-12) 

which  yields  17.7  cm  of  oil. 

The  water  table  was  maintained  5  cm  above  the  fine-coarse  interface.  Once  the 
water  table  was  established,  the  boundary  between  the  unsaturated  and  saturated  zones 
was  easily  distinguished.  In  order  to  exceed  the  maximum  predicted  thickness,  the  height 
of  oil  in  the  well  was  maintained  slightly  more  than  30  cm  above  the  interface  or  90  cm 
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above  the  bottom  of  the  flume.  The  LNAPL  was  injected  in  the  left  well  using  a  peristaltic 
pump.  Initially,  the  LNAPL  was  injected  at  90  ml/min,  then  the  pump  was  turned  on  and 
off  periodically  to  maintain  the  desired  oil  level  in  the  well.  After  injection,  the  mound 
height  was  allowed  to  decrease  until  the  LNAPL  lens  came  to  equilibrium.  The  maximum 
height  of  the  mound  measured  at  the  well  at  equilibrium  (110  hours)  was  17.7  cm,  which 
was  the  maximum  thickness  predicted.  Figures  4-6  a-h  show  the  computer  images  taken 
during  the  experiment.  The  computer  images  show  clearly  the  hysteresis  effect  at  point  A. 
As  the  LNAPL  shifts  from  an  LNAPL-water  drainage  curve,  to  an  LNAPL-water 
imbibition  curve,  the  bottom  of  the  lens  moves  upward.  The  bottom  of  the  mound  moved 
up,  and  pivoted  around  point  A  as  the  leading  edge  advanced  (point  B).  Figure  4-7  shows 
the  difference  in  oil  pressure  in  the  well  and  water  pressure  at  the  elevation  of  the  textural 
interface,  i.e.  at  the  advancing  edge  of  the  lens.  After  about  1 8  hours,  this  pressure 
difference  was  constant  at  the  entry  pressure  of  LNAPL-water  drainage  curve. 
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Figure  4-6.  Preliminary  experiment  showing  progression  of  lens  in  a  two-phase  system. 
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Figure  4-7,  Change  in  pressure  with  time. 

Three-fluid  phase  experiments  simulating  an  LNAPL  spill  were  conducted.  The 
objective  of  the  experiments  was  to  verify  the  assumptions  governing  the  infiltration  and 
equilibrium  position  and  thickness  of  the  lens.  Table  4-2  summarizes  the  results  of  the 
three-phase  experiments.  Notation  was  given  in  figure  4-5. 

The  flume  was  dry  packed  with  the  sieved  sand.  After  saturating  the  flume,  the 
water  levels  in  the  wells  were  lowered  to  44  cm  above  the  bottom  of  the  flume.  The  water 
was  allowed  to  drain  for  22  hours,  creating  a  capillary  fringe  height  of  21  cm.  A 
cumulative  volume  of  165  ml  of  fluid  set  1  was  injected  approximately  13  cm  above  the 
top  of  the  fringe.  The  injection  of  the  LNAPL  occurred  over  a  ten  minute  period.  After 
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the  injection  of  the  LNAPL  ceased,  the  water  table  height,  mound  height,  and  mound 
length  were  recorded  with  time. 

The  second  experiment  was  conducted  to  observe  the  effects  of  adding  less 
solution  than  experiment  1.  After  saturating  the  flume,  the  water  levels  in  the  wells  were 
lowered  to  18.5  cm  above  the  bottom  of  the  flume.  The  water  was  allowed  to  drain  for  22 
hours,  creating  a  capillary  fringe  height  of  21  cm.  A  cumulative  volume  of  100  ml  of 
Solution  1  was  injected  approximately  6  cm  above  the  top  of  the  fringe  (at  the  top  of  the 
sand).  The  injection  of  the  LNAPL  occurred  over  a  five  minute  period. 

After  equilibration  an  additional  50  ml  of  fluid  set  1  was  added  to  the  lens.  This  is 
referred  to  as  experiment  3.  It  was  assumed  that  during  the  infiltration  of  the  additional 
solution,  the  lens  would  spread  in  the  vertical  and  horizontal  directions  and  then 
rebound”  upward  to  its  final  position.  The  last  experiment  measured  the  final  thickness 
of  a  lens  using  fluid  set  2.  The  flume  was  dry  packed  with  the  F-14  sand.  The  top  of  the 
sand  was  at  a  height  of  90  cm  above  the  bottom  of  the  flume. 

Computer  images  taken  during  experiments  1  and  3  are  shown  in  figures  4-8  and 
4-9.  Graphs  of  lens  height  with  time  are  shown  in  figures  4-10,  4-1 1  and  4-12.  Measured 
and  predicted  values  are  given  in  Table  4-2.  The  final  measured  thicknesses  were  about 
15  cm.  Depending  on  whether  14  cm  or  15  cm  was  used  as  the  imbibition  entiy  pressure, 
the  predicted  thickness  was  17.8  cm  and  15.2  cm,  respectively.  The  value  for  AT  at 
equilibrium  was  approximately  2  cm  which  is  close  to  the  predicted  value  of  1.6  cm.  Thus, 
the  location  of  the  lens  with  respect  to  the  top  of  the  capillary  fringe  can  be  predicted  from 

the  interfacial  tensions  and  soil  properties.  Note  that  the  sum  of  AT  +  T  +  hw  should  be 
equal  to  21  cm. 
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Comparison  of  Measured  and  Predicted 
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Figure  4-12.  Mound  height  with  time  for  Experiment  3. 
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4.5  SUMMARY  AND  CONCLUSIONS 


A  conceptual  model  of  LNAPL  infiltration  and  spreading  at  the  capillary  fringe 

was  developed  and  verified  with  a  series  of  experiments  in  a  two-dimensional  flume. 

Conclusions  which  can  be  drawn  from  this  work  are 

•  the  top  of  the  lens  is  not  necessarily  at  the  top  of  the  capillary  fringe  but  its  location 
can  be  calculated  from  the  water-NAPL  and  NAPL-air  interfacial  tensions, 

•  capillary  pressures  at  each  point  on  the  water-NAPL  boundary  are  on  a  different 
hysteresis  loop  of  the  water-NAPL  capillary  pressure-saturation  curve, 

•  the  top  of  an  LNAPL  lens  is  flat  and  the  capillary  pressure  at  each  point  is  at  the  air 
entry  pressure  on  the  air-NAPL  imbibition  curve, 

•  the  capillary  pressure  at  the  advancing  edge  of  the  lens  is  the  entry  pressure  on  the 
water-NAPL  drainage  curve 

and 

•  the  maximum  possible  thickness  of  an  LNAPL  lens  is  the  difference  between  the 

capillary  pressure  entry  head  on  the  oil-water  drainage  curve  and  the  capillary  pressure 

entry  head  on  the  oil-water  imbibition  curve. 
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SECTION  III 


SECTION  III 


CHAPTER  1 

A  SEMI-ANALYTICAL  SOLUTION  FOR  ONE-DIMENSIONAL, 

TWO-FLUID  FLOW  IN  LAYERED  POROUS  MEDIA 

Abstract 

A  semi-analytical  solution  is  developed  for  simulating  transient,  one-dimensional, 
two-fluid  immiscible  flow  in  layered  media.  The  solution  uses  the  Runge-Kutta 
integration  method  to  solve  the  governing  flow  equations  and  an  iterative  shooting 
technique  to  insure  that  all  boundary  conditions  are  satisfied.  This  semi-analytical 
solution  extends  previous  two-fluid  solutions  to  vertical  flow,  layered  media,  and  finite- 
length  flow  regions.  Evaluation  of  two  initial-boundary-value  problems  suggests  that  the 
proposed  semi-analytical  solution  can  simulate  certain  physical  problems  in  two-fluid 
flow  without  experiencing  the  types  of  numerical  instabilities  commonly  associated  with 
traditional  finite  difference  and  finite  element  numerical  models.  The  primary 
disadvantages  of  the  solution  are  that  it  is  limited  to  one-dimensional  flow,  requires  the 
development  of  problem-specific  iteration  schemes  for  different  types  of  physical 
problems,  and  sometimes  does  not  converge  smoothly  within  a  time  step.  Continued 
development  of  this  semi-analytical  technique,  however,  will  likely  result  in  an  efficient 
method  for  solving  many  types  of  one-dimensional,  two-fluid  flow  problems. 

1.1  INTRODUCTION 

Mathematical  simulation  of  two-fluid  immiscible  flow  is  of  technical  interest  in 
many  fields  of  engineering  and  geology.  Practical  applications  of  two-fluid  flow  theory 
include  agricultural  infiltration  and  drainage  in  layered  soils,  migration  of  nonaqueous 
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phase  contaminants,  and  extraction  of  petroleum  products.  Simulation  of  two-fluid  flow 
in  layered  or  otherwise  heterogeneous  porous  media  has  commonly  been  difficult  using 
traditional  finite  difference  and  finite  element  models.  Due  to  the  highly  nonlinear 
governing  equations,  these  modeling  approaches  can  exhibit  numerical  instabilities  and 
convergence  problems,  particularly  when  sharp  saturation  fronts  move  across  material 
interfaces. 

A  semi-analytical  solution  has  been  developed  for  the  one-dimensional  Richards 
equation  which  does  not  appear  subject  to  the  types  of  stability  problems  that  can  affect 
traditional  numerical  models.  In  this  chapter,  the  general  solution  strategy  developed  for 
the  Richards  equation  is  extended  to  the  transient,  one-dimensional,  two-fluid  flow 
equations.  Initial-boundary-value  problems  in  two-fluid  flow  are  solved  using  the  Runge- 
Kutta  integration  technique  in  conjunction  with  the  shooting  method  to  satisfy  the 
specified  boundary  conditions.  This  approach  differs  from  finite  difference  and  finite 
element  models  in  the  following  ways: 

•  A  set  of  equations  is  solved  explicitly  and  sequentially,  whereas  finite  difference 
and  finite  element  formulations  generally  require  simultaneous  solution  of  a  set 
of  equations. 

•  The  spacing  between  nodes  is  varied  within  and  between  time  steps  based  on 
accuracy  criteria.  This  differs  from  traditional  models  which  rely  on  a  fixed-in¬ 
space  numerical  mesh. 

•  The  iteration  process  requires  convergence  of  only  one  or  two  boundary 
condition  values,  rather  than  iterations  being  performed  simultaneously  on 
values  associated  with  all  internal  nodes. 

Exact  integral  solutions  for  one-dimensional,  two-fluid  flow  have  been  developed 
by  McWhorter  and  Sunada  (1990).  These  solutions  are  limited  to  horizontal  flow 
problems  and  assume  a  semi-infinite,  homogenous  porous  medium.  The  semi- 
analytical  technique 
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described  in  this  chapter  extends  the  McWhorter-Sunada  solutions  to  vertical  flow, 
layered  systems,  and  finite-length  flow  regions.  The  resulting  transient  model  has  been 
successfully  applied  to  several  types  of  imbibition  and  drainage  problems  in 
homogeneous  and  layered  porous  media. 


1.2  THEORY 

Partial  differential  equations  for  one-dimensional  flow  of  two  immiscible, 
incompressible  fluids  have  been  derived  by  several  authors  (e.g.,  Buckley  and  Leverett, 
1942;  McWhorter,  1971;  Morel-Seytoux,  1973;  Fokas  and  Yortsos,  1982;  Corey,  1994). 

A  discussion  of  these  equations  and  associated  assumptions  is  provided  in  McWhorter 
and  Sunada  (1990). 

Referring  to  Figure  1-1 ,  the  flux  equation  for  two-fluid  flow  can  be  written  as: 

3S 

qw  =  qtf  -  sin(0)  E  -  D— (1-1) 

oX 

where  the  nomenclature  is  defined  in  Table  1.  In  this  form,  the  flux  equation  is  written 
using  effective  saturation  (Se)  as  the  dependent  variable.  Effective  saturation  is  a 
normalized  saturation  for  the  wetting  fluid  defined  by: 

w  ■  ^  d-2> 

Se  is  functionally  related  to  capillary  pressure  head  (hc),  which  is  defined  as  the 
difference  between  the  nonwetting  and  wetting  fluid  pressure  heads: 

hc  =  hnw  -  hw  (1-3) 
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Boundary  Condition(s) 


Figure  1-1.  One-Dimensional  Flow  Domain  for  the  Two-Fluid  Solution 


Table  1-1.  Nomenclature  Used  for  the  Two-Fluid  Solution. 

A 

C  Effective  saturation  capacity  [L'  ] 

D  Hydraulic  diffusivity  function  for  two-fluid  flow  [L2  T 1  ] 

E  Gravity-related  function  for  two-fluid  flow  [L  T"  ] 
f  Mobility  ratio  [dimensionless] 

g  Acceleration  of  gravity  [L  T*] 

h  Fluid  pressure  head  [L] 

hc  Capillary  pressure  head  [L] 

hd  Brooks-Corey  displacement  pressure  head  [L] 

Nonwetting-fluid  pressure  head  [L] 
hw  Wetting-fluid  pressure  head  [L] 

k  Intrinsic  permeability  of  porous  media  [L  ] 

Relative  permeability  of  nonwetting  fluid  [dimensionless] 

k^  Relative  permeability  of  wetting  fluid  [dimensionless] 

^nw  Nonwetting  fluid  flux  (Darcy  velocity)  [L  T^] 
qw  Wetting  fluid  flux  (Darcy  velocity)  [L  T-1] 
qt  Total  fluid  flux  (Darcy  velocity)  [L  T^] 

P  Fluid  pressure  [M  L-^  T2] 

Se  Effective  saturation  [dimensionless] 

Se‘  Effective  saturation  at  time  level  t  (end  of  last  time  step)  [dimensionless] 

S*+it  Effective  saturation  at  time  level  t+At  (end  of  current  time  step)  [dimensionless] 
S_  Volumetric  nonwetting-fluid  saturation  [dimensionless] 

S0  Initial  effective  saturation  at  t  =  0  [dimensionless] 

Sw  Volumetric  wetting-fluid  saturation  [dimensionless] 

Sr  Irreducible  volumetric  saturation  [dimensionless] 

t  Elapsed  time  [T] 

At  Duration  of  time  step  [T] 

x  Distance  from  boundary  [L] 

ax  Distance  between  two  adjacent  nodes  [L] 

0  Orientation  of  flow  system  (angle  counterclockwise  from  horizontal)  [degrees  or 
radians] 

<t>  Total  porosity  [dimensionless] 
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Table  1-1.  (continued) 


A  Brooks-Corey  pore  size  distribution  index  [dimensionless] 
p0  Reference  density  [M  L"3] 

Pnw  Density  of  nonwetting  fluid  [M  L'3] 
pw  Density  of  wetting  fluid  [M  L‘3] 

Mnw  Viscosity  of  nonwetting  fluid  [M  L"1  T-1] 

Mw  Viscosity  of  wetting  fluid  [M  L"**  T**] 

Sr  in  equation  (2)  is  the  irreducible  (minimum)  wetting  fluid  saturation  which  is 
approached  at  very  high  capillary  pressure  heads. 

In  the  flux  equation,  the  functions  f,  D,  and  E  are  defined,  respectively,  as: 


f(Se)  = 


+  kmw(S6)  Mwj-1 

krw(^e)  Mnw 


(1-4) 


D(Se)  = 


_  ( kP0g> 

(  Mnw  J 


krnwCge)  f(Se) 
C(SJ 


E(Se) 


nw 


w  J 


(  kPo9l 


V  ^nw 


kmw(Se)  f(Se) 


(1-5) 


(1-6) 


where: 


(1-7) 


and  and  krnw  are  the  relative  permeabilities  of  the  wetting  and  nonwetting  fluids, 
respectively.  In  equations  (1-5)  and  (1-6),  p0  is  an  arbitrary  reference  density  which 
relates  fluid  pressures  to  pressure  heads: 


\ 
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(1-8) 


Total  volume  flux  (qt)  is  defined  as  the  sum  of  the  wetting  and  nonwetting  fluid  fluxes: 


Clt  =  %  +  <ln 


If  the  porous  medium  is  rigid  and  contains  two  incompressible  fluids,  continuity 
relationships  lead  to: 


(1-10) 


which  requires  that  qt(x,t)  =  q,(t).  In  other  words,  at  any  given  time,  the  total  flux  is 
spatially  uniform  throughout  the  flow  system. 

In  terms  of  effective  saturation,  the  differential  form  of  the  wetting  fluid  continuity 
equation  in  the  x  direction  is: 


^  =  -<j)(1-Sr)^ 

ax  r  at 


(1-11) 


Taking  the  derivative  of  (1-1)  and  combining  with  (1-1 1)  leads  to  the  governing  flow 
equation  for  two-fluid  flow: 


a  Las, 


aV  +  sin(0) Us:  -q‘ 


df  Mas. 


dSJ  dx 


■  *P  -S.)3T  f1'12' 


Note  that  in  order  to  preserve  continuity  of  the  nonwetting  fluid,  the  two-fluid  model 
requires  that  differential  equations  (1-10)  and  (1-12)  be  solved  simultaneously. 
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1.3  BROOKS-COREY  RELATIONSHIPS 


In  equation  (1-12),  the  functional  relationships  between  capillary  pressure  head, 
effective  saturation,  and  relative  permeabilities  are  completely  arbitrary.  In  this  study, 
we  use  equations  developed  by  Brooks  and  Corey  (1964, 1966)  to  define  these 
constitutive  relationships.  Effective  saturation  is  related  to  capillary  pressure  head  as 


follows: 


8.00  = 


hcrx 

t~  I  for  h„  >  h„ 


(1-13) 


1  for  hc  *  hd 


where  X  is  the  pore  size  distribution  index  for  the  porous  medium  and  hd  is  the 
displacement  pressure  for  the  medium  and  fluids  of  interest.  Note  that  during  drainage, 
the  capillary  pressure  head  must  exceed  the  displacement  pressure  head  before 
desaturation  (Se  <  1 )  takes  place. 


Relative  permeabilities  for  the  wetting  and  nonwetting  fluids  are: 


USJ  -  S’'2'*  (1-14) 

-  (1  -S.)2  (l -S'™)  (1-15) 

Note  that  Corey  (1994)  presents  some  modifications  to  the  above  equation  for  kmw 
which  are  not  utilized  in  this  chapter.  The  function  C  is  given  by: 

C(S.)  -  -A  s.’*“  (1-16) 

hd 

All  other  functions  used  in  this  chapter  can  be  derived  by  algebraic  manipulation  and/or 
differentiation  of  the  above  equations. 
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1.4  NUMERICAL  APPROACH 

1.4.1  Approximation  of  the  Two-Fluid  Flow  Equation 

For  conditions  where  Se  -<  1 ,  we  approximate  the  time  derivative  in  (1- 
12)  using  a  finite  difference  approach: 


where  At  is  the  duration  of  the  current  time  step,  and  S0* ,  S‘+4t  are  the  effective 
saturations  at  x  occurring  at  the  beginning  and  end  of  the  time  step,  respectively.  Since 
time  (t)  is  no  longer  an  independent  variable,  we  express  the  left-hand  side  of  (1-17)  in 
terms  of  ordinary  derivatives.  Applying  the  chain  rule  to  the  left-hand  side  and  algebraic 
manipulation  leads  to: 


d2S, 


t+At 


dx2 


dD 

dS. 


/ 


Idx 

4>(1-Sr) 


df  1 
ds. 

-  sin(0) 

'  dE  ' 

dS‘+4t 

d 

dx 

At  D 


s';*'  + 


<t>(1-Sr) 

At  D 


Se‘(x) 


(1-18) 


Note  that  in  the  above  equation,  the  dependent  variable  Se  and  its  space  derivatives 
are  evaluated  at  time  level  t+At  (i.e.,  at  the  end  of  the  current  time  step),  which  implies 
a  fully  backward-in-time  numerical  scheme.  Furthermore,  the  value  S*(x)  is  known  a 
priori  from  the  results  of  the  last  time  step.  Adopting  the  notation,  Se  =  S*+it,  equation 
(1-18)  has  the  general  form: 


s"  =  F(Se,Se/,qt,x)  =  V(Se')2 


+  VS' 


+  VS. 


U 


(1-19) 
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where  f1  and  f3  are  functions  of  Se,  f2  is  a  function  of  Se  and  qt,  and  f4  is  a  function  of 
Se  and  x.  For  each  time  step,  this  second-order  ordinary  differential  equation  is 
integrated  to  determine  both  Se(x)  and  Se'(x)  using  the  fifth-order  Runge-Kutta 
algorithms  presented  in  Press  et  al.  (1992).  These  algorithms  contain  an  adaptive 
stepsize  control,  which  modifies  the  distance  between  nodes  (ax)  based  on  user- 
specified  accuracy  criteria.  Since  the  node  coordinates  change  each  time  the  solution 
is  evaluated,  values  of  S0*(x)  must  be  interpolated  from  results  of  the  previous  time 
step.  This  is  performed  using  the  cubic  spline  interpolation  algorithm  presented  in 
Press  etal.  (1992). 

Calculations  proceed  node-by-node  beginning  at  x  =  0.  Evaluation  at  a  node 
only  requires  information  determined  from  previous  nodes  within  the  sequence.  As  a 
result,  the  numerical  scheme  is  explicit  and  requires  only  the  sequential  solution  of 
independent  equations.  Once  the  values  Se(Xj)  and  Sfc)  are  obtained  for  a  node 
located  at  coordinate  xjt  the  capillary  pressure  head  at  that  node  is  computed  based  on 
equation  (1-13): 


W  =  hd  [Se(Xj)]_1/x 


(1-20) 


and  the  wetting  fluid  flux  at  the  node  is  computed  directly  from  equation  (1): 


qw(*i)  =  q,f[Se(x,)]  -  sin(0)  E[Se(xi)]  -  D^)]  Sfo)  (1-21) 


The  wetting  fluid  pressure  head  at  the  node  is  computed  using  a  finite  difference 
approximation  to  the  Darcy  equation: 


hw(xi)  =  hw(Xi-i) 


aXMw  I  qw(Xj-i)  A  qw(Xi)  1  AXpwsin(0) 
2kP0glkJSe(xi.1)]  kJSe(xi)]J  ~  p0 


(1-22) 
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where  ax  is  the  distance  between  nodes  Xj_^  and  x,.  Nonwetting  fluid  pressure  head  is 
computed  from  the  definition  of  capillary  pressure  in  equation  (1-3): 


h^)  =  he(Xj)  +  hw(Xj)  (1-23) 

By  equation  (1-2),  the  wetting  fluid  saturation  is  given  by: 

SJXj)  =  Sr  +  (1  -  Sr)  Se(Xj)  (1-24) 

and  based  on  continuity,  the  nonwetting  fluid  saturation  is: 


S^)  =  1  -  Sw(Xj) 


(1-25) 


1.4.2  Saturated  Conditions 

If  the  porous  medium  is  saturated  between  two  adjacent  nodes,  Se  =  1, 

=  1 ,  and  qw  has  the  same  value  at  both  nodes.  Equation  (1-22)  therefore  simplifies 
to: 


hJXj)  =  hjx,.,)  -  ax 


kpDg 


+ 


PwSin(e) 

Po 


(1-26) 


In  saturated  media,  the  nonwetting  fluid  does  not  exist  as  continuous  (mobile)  fluid  and 
thus,  its  pressure  head  does  not  have  physical  relevance  with  regard  to  fluid  flow. 


1.4.3  Material  Interfaces 

Referring  to  Figure  1-2,  at  the  interface  between  two  materials,  there  can 
be  a  discontinuity  in  effective  saturation,  but  pressure  heads  and  fluxes  must  be 
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Interface 


Flux  Continutiy 


Pressure  Continuity 


Saturation  Discontinuity  SeA  5*  SeB 


Figure  1-2.  Numerical  Formulation  Across  a  Material  Interface  for  the  Two-Fluid 
Solution. 
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continuous  across  the  interface.  Using  the  Brooks-Corey  relationships  and  taking 
advantage  of  capillary  pressure  head  continuity: 


S 


B 

e 


for  hcA  >  hdB 


1  for  hcA  <  hdB 


(1-27) 


where  the  superscript  A  indicates  the  capillary  pressure  head  computed  in  material  A  at 
the  interface,  the  result  of  previous  calculations,  and  B  indicates  capillary  properties 
associated  with  material  B.  Note  that  the  above  equation  is  applied  only  when 
evaluating  the  first  node  (boundary)  of  a  new  material  layer. 

If  the  magnitude  of  the  nonwetting  fluid  flux  (qnw)  is  greater  than  zero,  flux 
continuity  can  be  maintained  across  an  interface  only  if  the  effective  saturation  is  less 
than  unity  on  both  sides  of  the  interface.  When  the  Runge-Kutta  solution  is  extended 
across  an  interface,  there  can  exist  a  situation  where  the  magnitude  of  the  nonwetting 
fluid  flux  is  significant,  but  equation  (1-27)  computes  SeB  equal  to  one.  This  occurs 
when  the  capillary  pressure  in  material  A  is  less  than  the  displacement  pressure  for 
material  B.  In  this  case,  the  algorithm  artificially  increases  the  nonwetting  fluid  pressure 
on  the  B  side  of  the  interface  so  that  effective  saturation  (SeB)  is  reduced  to  0.999, 
which  is  equivalent  to  raising  the  capillary  pressure  to  just  slightly  above  the 
displacement  pressure  in  material  B.  This  procedure  allows  nonwetting  fluid  flux 
continuity  to  be  maintained,  but  creates  a  step-change  in  the  nonwetting  fluid  pressure 
head  and  capillary  pressure  head  across  the  interface.  Laboratory  experiments  have 
indicated  that  when  a  nonwetting  fluid  exits  the  outflow  face  of  a  porous  material,  there 
is  a  tendency  for  the  wetting  fluid  saturation  to  approach  unity  at  the  face  (Perkins, 
1957;  McWhorter,  1971;  Davis,  1990).  For  steady  flow  simulations,  Corey  (1994) 
recommends  that  capillary  pressure  head  at  the  outflow  face  be  specified  as 
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approximately  equal  to  the  displacement  pressure  head  of  the  porous  medium.  For 
material  interfaces  within  a  flow  region,  we  follow  this  recommendation  by  setting  the 
capillary  pressure  head  in  the  upstream  material  to  just  above  its  displacement  pressure 
head. 


1.4.4  Boundary  Conditions 

A  characteristic  of  the  numerical  formulation  is  that  the  solution  is 
constructed  explicitly,  beginning  at  x  =  0  and  then  proceeding  node-by-node  across  the 
flow  region.  To  begin  the  calculations,  an  input  value  of  At  and  boundary  conditions  at 
x  =  0  must  be  provided.  These  boundary  conditions  are  qw(0,t),  Se(0,t),  and  q,(t) .  Note 
that  effective  saturation  is  uniquely  defined  by  specifying  either  the  capillary  pressure 
head  or  the  pressure  heads  of  both  fluids.  These  values  can  be  used  in  place  of 
effective  saturation  as  a  boundary  condition.  Also,  for  incompressible  fluids,  the  total 
fluid  flux  at  time  t  is  a  constant  for  all  values  of  x.  When  using  the  Runge  Kutta  method, 
it  is  not  known  in  advance  what  values  of  qw,  qnw,  or  Se  will  be  computed  at  the 
opposite  (x  =  L)  boundary  of  the  system.  If  boundary  conditions  are  specified  at  x  =  L, 
an  iterative  procedure,  commonly  referred  to  as  the  shooting  method  (Press  etal., 

1992;  Acton,  1990;  Keller,  1968),  is  used  to  satisfy  those  conditions.  For  example, 
suppose  that  the  prescribed  boundary  conditions  for  the  problem  are  qt(t),  Se(0,t),  and 
Se(L,t)  as  shown  on  Figure  1-3.  In  this  case,  we  input  At,  qt(t),  the  prescribed  value  of 
Se(0,t),  and  a  trial  value  of  qw(0,t).  Then  qw(0,t)  is  varied  in  a  trial-and-error  manner 
until  the  computed  Se  at  the  opposite  end  of  the  flow  system  corresponds  to  the 
prescribed  x  =  L  boundary  condition.  Depending  on  the  problem  under  consideration, 
different  combinations  of  boundary  conditions  at  x  =  0  can  be  evaluated  using  the 
shooting  method  to  satisfy  a  prescribed  boundary  condition  at  x  =  L. 
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Prescribed  Boundary  Condition 
Trial  Value 
Computed  Value 
Trial  Solution 


True  (Converged)  Solution 


Figure  1-3.  Shooting  Method  Used  for  the  Two-Fluid  Solution. 
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1.5  EXAMPLE  SIMULATIONS 


1.5.1  Example  3A  -  Horizontal,  Two-Direction  Displacement  of  a 

Nonwetting  Fluid 

This  example  considers  the  displacement  of  a  nonwetting  fluid  by 
injection  of  wetting  fluid  at  one  end  of  a  horizontal,  finite-length  soil  column.  A 
description  of  this  initial-boundary-value  problem  is  shown  on  Figure  1-4  and  input 
parameters  are  summarized  in  Table  1-2.  Initially,  the  homogeneous  medium-grained 
sand  column  contains  two  fluids  under  no-flow  conditions.  The  initial  pressure  heads 
are  zero  for  the  wetting  fluid  and  32.32  cm  for  the  nonwetting  fluid.  At  the 
corresponding  capillary  pressure  head  of  32.32  cm,  the  initial  effective  saturation  (SQ)  is 
0.10  and  the  volumetric  wetting  fluid  saturation  is  0.19  throughout  the  column.  At  time 
zero,  the  wetting  fluid  pressure  head  at  x  =  0  increases  instantaneously  from  zero  to 
16.16  cm,  which  is  equivalent  to  instantaneously  raising  the  saturation  from  0.19  to  0.82 
at  the  left  side  of  the  column.  The  increase  in  wetting  fluid  pressure  and  saturation  at  x 
=  0  causes  the  wetting  fluid  to  imbibe  into  the  column  which  results  in  displacement  of 
the  nonwetting  fluid.  In  general,  the  displaced  nonwetting  fluid  will  exit  the  column  as 
counter-current  flow  at  x  =  0  and  unidirectional  flow  at  x  =  L;  in  other  words,  the 
nonwetting  fluid  displacement  process  is  two-directional.  The  relative  magnitudes  of 
nonwetting  fluid  outflows  at  each  end  of  the  column  depends  on  the  position  of  the 
saturation  front  over  time  and  the  pressure  heads  maintained  at  x  *  0  and  x  =  L. 

For  this  type  of  physical  problem,  a  two-step  iteration  procedure  is  used  for  each 
time  step.  First,  the  bisect  method  (Kreyszig,  1988)  is  used  to  iterate  on  qw(0,t),  based 
on  a  trial  input  value  of  the  total  fluid  flux  [qt(t)j.  Prior  to  the  saturation  front  reaching  the 
x  =  L  boundary,  these  iterations  are  performed  until  Se(x  *)  =  S0  and  Se'(x  *)  =  0  at 
some  coordinate  x*-;  L.  After  the  saturation  front  reaches  the  x  =  L  boundary,  the 
iterations  are  performed  until  Se(L)  =  S0 .  When  either  of  the  above  conditions  are  met, 
the  wetting  fluid  pressure  head  is  computed  at  x  =  L  and  compared  with  the  prescribed 
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^nw(®’p  32.32  cm  \  ^  /q 

hw(0,t)  =  16.16  cm  J  ^u,x' 


Initial  Condition  (t=0) 


Figure  1-4.  Example  3A  -  Horizontal,  Two-Direction  Displacement  of  a 
Nonwetting  Fluid  Solved  by  the  Two-Fluid  Solution. 
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Table  1-2. 


Parameters  Used  in  Example  Simulations  Based  on  the  Two-Fluid 
Solution. 


Pnw 

0.8 

g  cm*3 

Pw 

1 

g  cm*3 

(water) 

Po 

1 

gem*3 

(water) 

g 

980 

-2 

cm  s  * 

Pnw 

0.005 

-1  -1 
g  cm  1  s  1 

Mw 

0.01 

-1  -1 
g  cm  1  s  1 

e 

0 

degrees 

(Example  3A) 

90 

degrees 

(Example  3B) 

Coarse  Sand 

Medium  Sand 

k 

10‘6 

cm2 

ID’7  cm2 

X 

4 

3 

hd 

10 

cm 

15  cm 

sr 

0.05 

0.1 

0.35 

0.35 

boundary  condition  hw(L,t).  Based  on  the  difference  between  the  computed  and 
prescribed  wetting  fluid  pressures,  the  trial  value  of  qt(t)  is  adjusted  and  the  iteration 
process  repeated.  Convergence  for  the  time  step  is  achieved  when  the  computed 
wetting  fluid  pressure  head  at  x  =  L  corresponds  to  the  prescribed  boundary 
condition.  Using  the  secant  method  (Kreyszig,  1988),  convergence  is  usually 
achieved  with  three  to  six  adjustments  of  the  qt  value. 

Progression  of  the  saturation  front  into  the  column  is  shown  on  Figure  1-5.  Prior 
to  the  front  reaching  the  x  =  L  boundary,  the  saturation  profiles  have  a  characteristic 
shape  which  is  similar  to  that  of  nonwetting  fluid  displacements  in  a  semi-infinite  column 
(McWhorter  and  Sunada,  1990).  At  large  times,  after  the  front  has  reached  the  x  =  L 
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boundary,  the  saturation  profile  takes  on  the  characteristic  shape  for  steady-state  flow, 
in  this  case  with  the  nonwetting  fluid  flux  is  close  to  zero.  Figure  1-6  shows  saturation 
and  pressure  head  profiles  existing  at  an  elapsed  time  of  10,500  seconds.  Of  interest  is 
the  nonwetting  fluid  pressure  head  distribution  (Figure  1-6b).  Although  the  same 
nonwetting  fluid  pressure  head  is  prescribed  at  both  ends  of  the  column,  an  increased 
pressure  head  develops  inside  the  column.  This  creates  a  flow  divide  at  x  =  18.4  cm. 
Left  of  the  divide,  the  nonwetting  fluid  moves  toward  the  x  =  0  side  of  the  column 
(counter-current  to  the  wetting  fluid  flux)  and  right  of  the  divide,  nonwetting  fluid  flow  is 
toward  the  x  =  L  boundary  (same  direction  as  the  wetting  fluid  flux). 

The  magnitude  of  fluid  fluxes  at  the  boundaries  of  the  column  are  shown  on 
Figure  1-7.  At  early  times,  fluxes  occurring  at  x  =  L  are  relatively  small  and  wetting  fluid 
inflow  at  x  =  0  is  approximately  equal  to  the  nonwetting  fluid  outflow.  Therefore, 
essentially  full  counter-current  flow  occurs  at  this  early  stage  of  the  displacement 
process.  At  an  elapsed  time  of  about  7x10^  seconds,  the  nonwetting  fluid  outflow  at  x 
=  0  is  relatively  small  and  wetting  fluid  inflow  at  x  =  0  and  nonwetting  fluid  outflow  at  x  = 
L  are  approximately  equal.  This  indicates  that  displacement  process  is  nearly 
unidirectional  throughout  the  column. 

As  the  saturation  front  reaches  the  x  =  L  boundary,  there  is  a  rapid  increase  in 
the  wetting  fluid  outflow  and  at  large  times,  a  steady-state  condition  is  reached  where 
the  wetting  fluid  inflow  at  x  =  0  is  approximately  equal  to  its  outflow  at  x  =  L.  Because 
the  nonwetting  fluid  pressure  heads  are  the  same  at  both  ends  of  the  column,  the 
nonwetting  fluid  flux  should  tend  towards  zero  throughout  the  column  as  steady-state  is 
approached.  On  Figure  1-7,  this  is  clearly  shown  for  the  nonwetting  fluid  outflow  at  x  = 

L.  At  large  time,  the  nonwetting  fluid  outflow  at  x  =  0  is  shown  to  reach  a  steady-state 
flux  of  about  two  orders  of  magnitude  lower  than  the  wetting  fluid  fluxes.  This  computed 
nonzero  outflow  rate  appears  to  be  a  numerical  artifact  related  to  the  overall  accuracy  of 
the  semi-analytical  solution. 
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Wetting  Phase  Saturation: 


* 
c n 


Figure  1-5.  Example  3A  -  Saturation  Profiles. 
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hw(cm  water)  hnw  (cm  water) 


1.0 


a.  Saturation 


Figure  1-6.  Example  3A  -  Saturation  and  Pressure  Head  Profiles, 
t  =  10,500  seconds 
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-  Wetting  Fluid  Inflow  at  x  =  0: 

- Nonwetting  Fluid  Outflow  at  x  =  0: 

- Wetting  Fluid  Outflow  at  x  =  L: 

- Nonwetting  Fluid  Outflow  at  x  =  L: 


q„(L.t) 

<w(m) 


Figure  1-7.  Example  3A  -  Wetting  and  Nonwetting  Fluid  Fluxes. 
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1.5.2  Example  3B  -  Vertical  Migration  of  LNAPL  in  a  Layered  Soil 

This  example  considers  the  downward  migration  of  a  light  nonaqueous 
phase  liquid  (LNAPL),  through  a  layered  soil,  in  response  to  a  lowering  of  the  water 
table.  This  constitutes  a  drainage  problem  where  the  wetting  fluid  (water)  is  displaced 
by  the  nonwetting  fluid  (LNAPL).  The  physical  problem  for  the  simulation  is  shown  on 
Figure  1-8  and  input  parameters  are  summarized  in  Table  1-2.  The  flow  domain 
consists  of  a  three-layer  system  with  a  zone  of  coarse  sand  situated  between  two  layers 
of  medium  sand.  Initially,  the  system  is  fully  saturated  with  the  water-table  located  at 
the  top  of  the  soil  profile.  For  static  conditions,  this  results  in  a  water  pressure  head  of 
100  cm  at  the  base  of  the  profile.  At  time  zero,  the  water  pressure  head  at  the  base  of 
the  profile  is  instantaneously  lowered  from  100  to  67  cm,  and  a  shallow  ponding 
condition  for  the  LNAPL  is  imposed  at  the  top  of  the  profile.  Since  the  fluids  and  porous 
media  are  assumed  to  be  incompressible,  the  water  pressure  head  at  the  top  of  the 
profile  decreases  at  time  zero  and  the  associated  capillary  pressure  head  becomes 
greater  than  displacement  pressure  head.  As  a  consequence,  water  drains  at  the  top  of 
the  soil  profile  and  is  displaced  by  downward  moving  LNAPL.  The  displaced  water  exits 
the  system  as  outflow  at  the  bottom  of  the  profile. 

For  this  type  simulation,  a  two-step  iteration  procedure  is  used  for  each  time 
step.  We  define  a  desaturation  point  (x  *)  such  that  the  effective  saturation  (Se)  is 
equal  to  1  for  a\\  x  <  x  *  and  is  less  than  1  for  x  >  x  *.  The  first  step  is  to  input  a  trial 
value  of  qw(0,t)  and  then  iterate  on  the  x*  value.  For  x  -<  x*.  the  system  is  fully 
saturated,  qw  =  qt  =  constant,  and  the  water  pressure  head  (hw)  distribution  is 
computed  using  the  saturated  form  of  Darcy's  law.  At  x  =  x*,  the  following  conditions 
exist:  Se  =  1,  Se'  -(-«),  and  qt  =  qw(0,t).  These  conditions  are  used  as  the  starting 
point  for  the  Runge-Kutta  method  to  compute  the  saturation  distribution  within  the 
x  >  x  *  region.  Using  the  bisect  method,  iterations  are  performed  on  the  x  *  value  until 
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Medium  Sand 
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gogj  Coarse  Sand 


Boundary  Conditions  (t  >  0) 


hnw(Ut)  =  0 

qw(L,t)  =  0 


Initial  Condition  (t=0) 
Se(x,0)  =  1 


Boundary  Conditions  (t  >  0) 
hw(0,t)  =  67  cm 


Figure  1-8.  Example  3B  -  Downward  Migration  of  LNAPL  in  a  Layered  Soil 
Solved  by  the  Two-Fluid  Solution. 


q^L.t)  is  computed  to  be  approximately  equal  to  zero,  thus  satisfying  the  prescribed  no¬ 
water-flow  boundary  condition  at  the  top  of  the  profile.  The  next  step  is  to  compute  the 
nonwetting  fluid  pressure  head  at  the  top  of  the  profile  [hnw(L,t)].  To  satisfy  the  shallow 
ponding  condition  for  LNAPL,  this  value  should  be  zero.  If  this  value  is  significantly 
different  than  zero,  the  qw(0,t)  input  value  is  adjusted  and  iterations  are  once  again 
performed  on  x  *.  Using  the  secant  method,  the  qw(0,t)  input  value  typically  requires 
three  to  six  adjustments  to  satisfy  the  prescribed  hnw(L,t)  boundary  condition.  Under 
certain  conditions,  the  desaturation  point  may  exist  at  a  material  interface  and  the 
effective  saturation  at  the  interface  falls  below  one.  In  this  case,  x*  is  fixed  at  the 
interface  coordinate  and  iterations  are  performed  on  Se  rather  than  x*. 

The  progression  of  saturation  profiles  over  time  is  shown  on  Figure  1-9.  At  early 
times  (Figure  1-9a),  LNAPL  has  migrated  downward  only  in  the  upper  medium  sand 
layer  and  the  desaturation  point  (x  *)  is  located  at  x  =  85  cm.  At  a  time  of  4,040 
seconds,  the  saturation  front  reaches  the  first  interface  and  exhibits  a  saturation 
discontinuity  as  shown  on  Figure  1-9b.  The  wetting  fluid  saturation  is  nearly  equal  to 
one  just  above  the  interface  and  is  significantly  less  than  one  below  the  interface.  As 
time  progresses  (Figures  1-9c  and  1-9d),  the  saturation  profile  is  relatively  stable  in  the 
upper  layer  and  most  desaturation  occurs  in  the  middle  coarse  sand  layer.  Note  that  at 
12,040  seconds  (Figure  1-9d),  the  bottom  of  the  upper  layer  is  just  beginning  to 
desaturate,  but  the  saturation  discontinuity  across  the  interface  remains.  Also,  a 
saturation  discontinuity  has  developed  across  the  interface  between  the  middle  and 
lower  layers.  At  14,190  seconds  (Figure  1-9e),  the  middle  coarse  sand  layer  is  mostly 
desaturated,  the  front  has  moved  past  the  second  interface,  and  continued  desaturation 
occurs  at  the  bottom  of  the  upper  layer.  With  increased  time  (Figures  1-9f  and  1-9g), 
slow  desaturation  takes  place  primarily  in  the  upper  and  lower  layers.  At  infinite  time 
(Figure  1-9h),  the  desaturation  process  ceases  and  a  new  static  condition  is 
established.  At  this  stage,  there  is  no  further  LNAPL  migration  and  all  fluxes  are  zero. 
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The  variation  in  the  total  fluid  flux  (qt)  over  time  is  shown  on  Figure  1-10.  For 
this  problem,  the  magnitude  of  the  total  fluid  flux  is  equal  to  the  LNAPL  inflow  rate  at  the 
top  of  the  soil  profile  and  the  water  outflow  rate  at  the  bottom  of  the  profile.  After  the 
saturation  front  reaches  the  first  interface  (about  4,000  seconds)  there  is  a  temporary 
increase  in  the  total  fluid  flux,  which  is  presumably  caused  by  the  relatively  rapid  initial 
desaturation  of  the  middle  coarse  sand  layer.  Between  about  4,000  and  12,000 
seconds,  the  flux  decreases  gradually  as  the  middle  layer  continues  to  desaturate.  At 
about  12,000  seconds,  the  front  interacts  with  the  second  interface  and  there  is  a  sharp 
decrease  in  the  total  fluid  flux.  At  large  times,  the  middle  layer  is  almost  completely 
desaturated  and  very  slow  desaturation  takes  place  in  the  upper  and  lower  layers.  As  a 
consequence,  the  flux  decreases  to  small  values. 


1.6  DISCUSSION 

The  Runge-Kutta  and  shooting  methods  employed  in  this  solution  do  not  exhibit 
the  types  of  stability  problems  commonly  associated  with  traditional  finite  difference  and 
finite  element  numerical  models.  We  have  not  experienced  any  stability  problems  when 
a  saturation  front  crosses  a  material  interface,  even  if  large  discontinuities  exist  in  both 
Se  and  Se'  at  the  interface.  The  proposed  solution  method  seems  capable  of  efficiently 
solving  certain  types  of  two-fluid  flow  problems  that  have  been  difficult  to  simulate  with 
traditional  models. 

There  are  occasions,  however,  when  convergence  within  a  time  step  is  not 
easily  achieved  using  this  solution  technique.  Adjustment  of  accuracy  parameters 
and/or  the  time  step  size  nearly  always  achieves  convergence,  but  this  can  be 
something  of  a  trial-and-error  process.  We  have  intentionally  designed  the  current 
computer  code  with  minimum  automation  capabilities  in  selecting  or  adjusting  these 
numerical  parameters.  It  is  our  expectation  that  alternative  solution  methods,  improved 
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Total  Flux  Rate:  |qt|  (x  1CT3  cm/sec) 


Figure  1-10.  Example  3B  -  Total  Volume  Flux. 


111-28 


programming,  and  increased  numerical  precision,  will  greatly  enhance  the  convergence 
characteristics  of  this  solution  method. 

One  significant  limitation  of  this  solution  is  that  problem-specific  iteration 
schemes  must  be  designed  for  different  types  of  physical  problems,  and  sometimes  the 
iteration  scheme  itself  must  change  during  the  course  of  a  simulation.  Development  of 
efficient  iteration  strategies  can  be  time  consuming  and  usually  requires  some  physical 
understanding  and  intuition  of  the  physical  problem  under  consideration. 


1.7  CONCLUSIONS 

The  semi-analytical  solution  presented  in  this  chapter  appears  to  be  a  promising 
technique  for  solving  certain  types  of  one-dimensional  problems  in  two-fluid  flow.  Our 
goal  in  this  chapter  has  been  to  present  the  basic  solution  strategy  and  show  that  it  is 
capable  of  solving  problems  of  practical  interest  in  two-fluid  flow.  The  method  is 
particularly  well  suited  for  simulating  flow  in  layered  media  and  does  not  exhibit  the 
types  of  stability  problems  that  can  occur  in  traditional  finite  difference  and  finite  element 
numerical  models.  The  primary  disadvantages  of  the  proposed  solution  are  that  it  is 
limited  to  one-dimensional  flow,  requires  the  development  of  problem-specific  iteration 
schemes  for  different  types  of  physical  problems,  and  sometimes  does  not  converge 
smoothly  within  a  time  step.  Continued  development  of  this  semi-analytical  technique, 
however,  will  likely  result  in  an  efficient  method  for  solving  many  types  of  one¬ 
dimensional,  two-fluid  flow  problems. 
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CHAPTER  2 


COMPARISON  OF  THE  RICHARDS  EQUATION  AND  TWO-FLUID  FLOW 
EQUATIONS  FOR  DRAINAGE  WITH  LIMITED  AIR  ACCESS 


Abstract 

Analysis  of  drainage  in  air-water  systems  has  traditionally  been  based  on 
solutions  to  the  Richards  equation.  Utilizing  analytical  models  for  two-fluid  flow  with 
zero  and  finite  air  viscosities,  the  assumptions  of  the  Richards  model  are  critically 
evaluated  for  the  case  where  air  access  into  the  porous  medium  is  restricted.  It  is 
shown  that  for  certain  types  of  drainage  problems,  the  Richards  model  does  not  provide 
for  accurate  simulation  of  drainage  behavior,  particularly  in  layered  soils  with  real  air 
viscosity.  Furthermore,  there  is  a  class  of  drainage  problems  for  which  the  Richards 
equation  does  not  allow  a  theoretical  solution  to  be  obtained.  In  these  cases,  rigorous 
simulation  of  the  drainage  process  can  only  be  performed  by  considering  the  two-fluid 
(air-water)  dynamics  of  the  system.  Because  the  Richards  equation  assumes  that  no  air 
pressure  gradients  exist  within  the  system,  it  could  be  interpreted  to  contain  an  implicit 
assumption  of  zero  air  viscosity.  This  conceptualization  is  inaccurate  and  potentially 
misleading,  because  the  Richards  equation  is  not  recovered  from  the  two-fluid  flow 
equations  in  the  limit  of  zero  air  viscosity.  It  is  more  appropriate  to  view  the  Richards 
equation  as  an  unlimited  air  access  model,  where  air  is  supplied  to  the  system  at  zero 
pressure  head. 


2.1  INTRODUCTION 

Mathematical  simulation  of  water  drainage,  under  partially  saturated  conditions, 
is  of  practical  interest  in  agricultural  and  geotechnical  engineering.  During  the  drainage 
process,  any  decrease  in  water  saturation  must  be  accompanied  by  movement  of  air 
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into  the  system  to  occupy  pore  space  previously  containing  water.  Drainage  analyses 
using  analytical  equations  or  numerical  models  have  traditionally  been  based  on 
solutions  to  the  Richards  equation  (Richards,  1931;  Jury  etal.,  1991).  Two  basic 
assumptions  in  the  Richards  equation  are  (1)  the  air  pressure  is  spatially  uniform  and 
(2)  air  is  freely  accessible  to  all  points  within  the  flow  system  regardless  of  the 
saturation  distribution.  Some  investigators  (Corey,  1994;  Dullien,  1992)  have 
considered  the  Richards  equation  to  be  applicable  to  systems  in  which  air  offers  no 
resistance  to  flow,  which  could  be  interpreted  to  imply  that  the  air  viscosity  is  effectively 
equal  to  zero.  At  face  value,  this  interpretation  may  seem  reasonable  because  air 
viscosity  is  only  about  2  percent  of  the  viscosity  of  water.  Furthermore,  at  zero 
viscosity,  air  movement  can  take  place  with  no  pressure  gradient,  which  is  consistent 
with  the  uniform  air  pressure  assumption. 

Technical  interest  in  the  simultaneous  flow  of  immiscible  fluids  has  led  to  the 
development  of  analytical  solutions  and  numerical  models  which  can  be  used  to  analyze 
two-fluid  (air-water)  flow  in  porous  media.  These  mathematical  models  can  explicitly 
account  for  viscous  air  flow  and  are  not  limited  by  the  uniform  air  pressure  assumption 
used  in  the  Richards  equation.  In  addition,  simulations  can  be  performed  that  consider 
limitations  on  the  access  of  air  into  the  medium.  With  the  ability  to  simulate  two-fluid 
flow  in  homogeneous  and  layered  systems,  it  is  of  interest  to  evaluate  the  applicability  of 
the  Richards  equation  for  simulating  practical  problems  of  drainage  in  air-water 
systems. 

Our  treatment  of  drainage  gives  full  consideration  to  the  simultaneous  bulk  flow 
of  air  and  water.  The  potential  effects  of  air  flow  on  the  drainage  of  water  have  not 
been  systematically  investigated.  Most  studies  of  simultaneous  air-water  flow  have 
focused  on  the  infiltration  problem  (McWhorter,  1971;  Morel-Seytoux,  1969,  1973).  An 
approximate  analysis  of  column  drainage  as  a  two-fluid  problem  is  presented  in  Morel- 
Seytoux  (1975). 
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In  this  chapter,  the  assumptions  of  the  Richards  equation  are  critically  evaluated 
considering  two-fluid  flow  theory.  Attention  is  focused  on  details  of  the  flux  equations 
associated  with  the  Richards  and  two-fluid  flow  models.  In  addition,  we  develop  and 
compare  a  third  analysis  approach,  termed  the  zero-viscosity  model,  which  is  derived 
from  the  two-fluid  equations  in  the  limit  as  the  air  viscosity  goes  to  zero.  Of  particular 
interest  is  the  applicability  of  the  various  models  to  drainage  problems  in  layered  soils. 

To  evaluate  initial-boundary-value  problems  for  both  the  Richards  equation  and  the 
two-fluid  flow  equations,  use  is  made  of  the  semi-analytical  solutions  developed  in 
Chapters  2  and  3  of  this  dissertation. 

2.2  FAILURE  OF  THE  RICHARDS  EQUATION  TO  SIMULATE  A  SIMPLE 

DRAINAGE  PROBLEM 

This  evaluation  has  been  prompted  by  the  fact  that  the  Richards  formulation  is 
apparently  unable  to  provide  a  theoretical  solution  to  certain  types  of  air-water  drainage 
problems.  An  example  of  such  a  problem  is  presented  in  this  section.  Referring  to 
Figure  2-1,  consider  vertical  drainage  in  a  laboratory  soil  column  having  impermeable 
walls.  Initially,  the  column  is  fully  saturated  and  hydrostatic,  with  the  phreatic  (zero 
water  pressure)  surface  located  at  the  top  of  the  soil.  At  time  zero,  the  valve  is  opened 
and  water  discharges  to  the  atmosphere  at  the  end  of  a  drainage  pipe.  Laboratory 
experience  indicates  that  under  these  conditions,  air  will  spontaneously  invade  the  top 
of  the  column  and  displace  the  water  which  drains  out  the  bottom.  With  time,  a 
desaturation  front  will  migrate  downward  from  the  top  of  the  column.  Because  the  sides 
and  base  of  the  column  are  sealed,  and  the  lower  portion  of  the  soil  is  saturated,  there 
is  no  air  access  to  the  bottom  of  the  column.  Thus,  no  desaturation  can  occur  at  the 
base  of  the  soil  until  such  time  that  the  desaturation  front  has  migrated  downward  to  this 
level. 
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Shown  on  Figure  2-1  is  a  saturation  profile  at  some  time  after  initiation  of  the 
drainage  process.  We  define  an  incipient  desaturation  point,  which  is  located  at 

coordinate  x  .  The  column  is  fully  saturated  below  the  desaturation  point  and  partially 
saturated  above  the  point. 

In  the  region  above  and  including  the  desaturation  point,  the  Richards  flux 
equation,  presented  in  Chapter  2,  can  be  written  for  vertical  flow  as: 


l<U  = K kJSJ 


1  _  i  as; 

CJSJ  dx  t 


for  Sw  *  1 


(2-1) 


where  |qw|  is  the  magnitude  of  the  (downward)  water  flux,  K  is  the  saturated  hydraulic 

conductivity  of  the  porous  medium,  is  the  relative  water  permeability,  and  saturation 
capacity  is  defined  as: 


OJ.SJ  -  ^ 

dhc 


(2-2) 


At  the  desaturation  point  (Sw  =  1),  the  relative  water  permeability  is  equal  to  unity  and 
equation  (2)  becomes: 


K*\  -  K 


1  _  1  ggw 

CJ1)  ax  ; 


at  x  =  x* 


(2-3) 


Saturation  capacity  always  has  a  negative  value  and  the  saturation  profile  on  Figure  2-1 
shows  that  at  x  *,  the  derivative  of  Sw  must  be  either  zero  or  a  negative  value.  With 
these  considerations,  it  must  be  concluded  that  for  the  Richards  model: 


l^lw*lmax  ~  ^  Richards  Model 


(2-4) 
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Figure  2-1.  Vertical  Drainage  in  a  Homegeneous  Soil  Column. 
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This  equation  indicates  that  at  the  desaturation  point,  the  maximum  water  flux  that  can 
be  computed  using  the  Richards  flux  equation  is  equal  to  the  saturated  hydraulic 
conductivity  of  the  material  containing  the  desaturation  point. 

Capillary  pressure  head  (hc)  is  defined  as  the  difference  between  the  air  and 
water  pressure  heads.  At  the  desaturation  point,  the  water  saturation  is  equal  to  unity 
and  the  Brooks-Corey  (1964, 1966)  saturation  relationships  require  that  the  capillary 
pressure  head  be  equal  to  the  displacement  pressure  head  (hd).  Using  the  traditional 
Richards  assumption  that  the  air  pressure  head  is  everywhere  equal  to  zero,  the 
capillary  pressure  head  becomes  equal  to  the  negative  of  the  water  pressure  head  (hw). 
The  above  considerations  lead  to: 


hw  =  -hd  at  x  =  x* 


(2-5) 


Because  the  column  is  saturated  for  x  *  x* ,  the  water  flux  at  the  desaturation  point  can 
be  computed  alternatively  using  the  saturated  form  of  Darcy's  law.  Considering  that  the 
water  pressure  head  is  equal  to  -hd  at  x  =  x*  and  -B  at  x  =  0,  the  water  flux  in  the 
saturated  portion  of  the  column,  including  the  desaturation  point,  is: 


|qw*l  =  K 


B~hd 

x* 


+  1 


Darcy’s  Law 


(2-6) 


The  dimension  B  is  equal  to  the  magnitude  of  the  water  pressure  (suction)  head  at  the 
base  of  the  soil.  It  is  observed  from  the  above  equation  that  if  the  dimension  B  is 
greater  than  the  value  of  the  displacement  pressure  head  (hd),  Darcy's  law  gives  a 
water  flux  at  x  =  x*  which  is  larger  than  the  value  that  can  be  theoretically  computed  by 
the  Richards  flux  equation  at  the  same  point.  From  a  physical  standpoint,  the  elevation 
of  the  discharge  pipe  is  completely  arbitrary.  Thus,  for  the  case  where  B  >  hd,  the 
physical  problem  shown  on  Figure  2r1  cannot  be  solved  using  the  Richards  formulation. 
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For  this  situation,  the  Richards  equation  simply  does  not  apply  because  it  is  not  capable 
of  maintaining  flux  continuity  at  the  desaturation  point  when  B  >  hd .  Under  these 
conditions,  the  magnitude  of  the  water  flux  at  the  desaturation  point  exceeds  the 
saturated  hydraulic  conductivity.  If  the  Richards  flux  equation  (3)  is  solved  for  the 
saturation  derivative  (dS Jdx),  it  will  predict  that  the  derivative  is  positive  and  greater 
than  zero  at  the  desaturation  point.  This  leads  to  a  contradiction  because  it  would 
erroneously  indicate  that  there  is  no  desaturation  directly  above  the  desaturation  point. 


2.3  GOVERNING  EQUATIONS  FOR  ONE-DIMENSIONAL  FLOW 

To  systematically  evaluate  the  above  failure  of  the  Richards  equation,  we 
present  three  mathematical  models  which  can  be  used  to  evaluate  air-water  drainage. 
The  governing  equations  for  these  models  are  discussed  below. 

2.3.1  General  Relationships 

The  mathematical  nomenclature  used  in  this  chapter  is  summarized  in 
Table  2-1.  We  write  the  governing  differential  equations  using  effective  saturation  (Se) 
as  the  dependent  variable.  Effective  saturation  is  a  normalized  water  saturation  defined 
by: 

s.w  ■  <2-7> 

where  Sw  is  the  volumetric  water  saturation.  Se  is  functionally  related  to  capillary 
pressure  head  (hc),  which  is  defined  as  the  difference  between  the  air  and  water 
pressure  heads: 


(2-8) 
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Table  2-1.  Nomenclature  Used  in  Chapter  2. 


C  Effective  saturation  capacity  [L-^] 

Cw  Water  saturation  capacity  [L"^J 

D  Hydraulic  diffusivity  function  used  in  two-fluid  model  [L2  T1] 

Dr  Hydraulic  diffusivity  function  used  in  the  Richards  model  [L2  T^] 

Dv  Hydraulic  diffusivity  function  used  in  zero-viscosity  model  [L2  T1] 

E  Gravity-related  function  used  in  two-fluid  model  [L  T"^] 

Er  Gravity-related  function  used  in  the  Richards  model  [L  r1] 

Ev  Gravity-related  function  used  in  zero-viscosity  model  [L  T"1] 
f  Mobility  function  used  in  two-fluid  model  [dimensionless] 
fv  Mobility  function  used  in  zero-viscosity  model  [dimensionless] 

g  Acceleration  of  gravity  [L  T2] 

h  Fluid  pressure  head  [L] 

ha  Air  pressure  head  [L] 

hc  Capillary  pressure  head  [L] 

hd  Brooks-Corey  displacement  pressure  head  [L] 

1^,  Water  pressure  head  [L] 

K  Saturated  hydraulic  conductivity  [L  T^] 

k  Intrinsic  permeability  of  porous  media  [L2] 

k,,  Relative  air  permeability  [dimensionless] 

k^  Relative  water  permeability  [dimensionless] 

qa  Air  flux  (Darcy  velocity)  [L  T1] 

q,  Total  volume  flux  (Darcy  velocity)  [L  T^] 

qw  Water  flux  (Darcy  velocity)  [L  T^] 

q*  Water  flux  at  desaturation  point  [L  T"1] 

Sa  Volumetric  air  saturation  [dimensionless] 

Se  Effective  saturation  [dimensionless] 

Sr  Irreducible  volumetric  saturation  [dimensionless] 

Sw  Volumetric  water  saturation  [dimensionless] 

t  Elapsed  time  [T] 

x  Distance  coordinate  [L] 

x  *  Coordinate  at  desaturation  point  [L] 
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Table  2-1.  (continued) 


0  Total  porosity  [dimensionless] 

A  Brooks-Corey  pore  size  distribution  index  [dimensionless] 

pw  Water  density  [M  L'3] 

pa  Air  viscosity  [M  T^] 

Water  viscosity  [ML  T  ] 


In  this  chapter,  all  pressure  heads  are  expressed  in  terms  of  the  equivalent  height  of 
water  in  a  manometer.  Sr  is  the  irreducible  (minimum)  water  saturation  which  is 
approached  at  very  high  capillary  pressure  heads.  Since  the  system  contains  only  two 
fluids,  the  volumetric  air  saturation  given  by: 

Sa  =  1  -  Sw  (2-9) 

In  terms  of  effective  saturation,  the  differential  form  of  the  water  continuity 
equation  in  the  x  direction  is: 


^  -  -<K1-S,)^l  (2-10) 

d  X  (7 1 

2.3.2  Brooks-Corey  Relationships 

The  functional  relationships  between  capillary  pressure  head,  effective 
saturation,  and  relative  permeabilities  are  completely  arbitrary.  In  this  study,  we  use 
equations  developed  by  Brooks  and  Corey  (1964,  1966)  to  define  these  constitutive 
relationships.  Effective  saturation  is  related  to  capillary  pressure  head  as  follows: 
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MO  - 


(2-11) 


[  1  for  hc  =s  hd 

where  A  is  the  pore  size  distribution  index  and  hd  is  the  air-water  displacement  pressure 
for  the  medium  of  interest.  Note  that  during  drainage,  the  capillary  pressure  head  must 
exceed  the  displacement  pressure  head  before  desaturation  (S#-<  1 )  takes  place. 
Relative  permeabilities  for  water  and  air  are  defined,  respectively,  as: 


USJ  -  S.9,M  (2-12) 

USJ  -  (1-S.)2(l-Ss',w)  (2-13) 

Note  that  in  Corey  (1994),  certain  modifications  are  presented  for  equation  (13)  which 
are  not  considered  in  this  chapter.  Effective  water  saturation  capacity  is  given  by: 

c«j  -  g  -  -A  sj*“  (2-14, 

Note  that  the  value  of  C(Se)  is  always  negative.  All  other  functions  used  in  this  chapter 
can  be  derived  by  algebraic  manipulation  and/or  differentiation  of  the  above  equations. 


2.3.3  Two-Fluid  Flow  Model 

Referring  to  Figure  2-2a,  the  flux  equation  for  vertical  two-fluid  flow  can 

be  written  as: 


%  ~  ~  ^ 


-  D 


as, 

ax 


where  the  functions  f,  D,  and  E  are  defined,  respectively,  as: 


(2-1 5a) 
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f(S.)  = 


M 


(2-1 5b) 


1  +  kra(Se)  Mw 

kJSe)  Ma 


D(S ) 


(  kPwQ'l 


l  M.  j 


US.)  f(S.) 


E(Se)  = 


(  kPw9^ 


l  J 


C(S.) 


US.)  f(Se) 


(2-1 5c) 

(2-1 5d) 


Taking  the  derivative  of  (2-15)  and  combining  with  (2-10)  leads  to  the  first  partial 
differential  equation  (PDE)  for  two-fluid  flow: 


a  1 

'  as.l 

D — - 

+ 

dE  1 

-  Pt 

[—11 

dx\ 

k  j 

ldSej 

[dS.jj 

as. 

ax 


<D(i  -sr) 


at 


(2-16) 


Total  volume  flux  (qt)  is  defined  as  the  sum  of  the  water  and  air  fluxes: 


Pt  -  %  +  Pa 


(2-17) 


If  the  porous  medium  is  rigid  and  contains  two  incompressible  fluids,  continuity 
relationships  lead  to  the  second  PDE  for  two-fluid  flow: 


dPt 

ax 


=  o 


(2-18) 


which  requires  that  q,(x,t)  =  qt(t) .  In  other  words,  at  any  given  time,  the  total  flux  is 
spatially  uniform  throughout  the  flow  system.  For  many  physical  problems  of  air-water 
flow,  the  assumption  of  an  incompressible  air  phase  is  not  reasonable.  For  drainage 
problems,  however,  we  consider  this  assumption  to  be  acceptable  as  a  first-order 
approximation.  Note  that  the  two-fluid  model  requires  that  the  two  PDE’s  in  (16)  and 
(18)  be  solved  simultaneously. 
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a.  Two— Fluid 
Model 


b.  Richards 
Model 


Figure  2-2.  One-Dimensional  Drainage  Flow  Models. 


2.3.4  Richards  Model 


Referring  to  Figure  2-2b,  the  Richards  flux  equation  for  partially  saturated 
water  flow  can  be  written  as: 


(2-1 9a) 


where  the  functions  DR  and  ER  are  defined,  respectively  as: 


DR(Se)  =  - 


kpw9  krw(Se) 
Mw  C(Se) 


(2-1 9b) 


Er(S6)  =  ^  USe) 
Hw 


(2-1 9c) 


Taking  the  derivative  of  (2-19)  and  combining  with  (2-10)  leads  to  the  following  PDE: 


d 

kas*l 

JdE‘) 

dX 

ax  J 

ldSeJ 

as, 

dx 


(2-20) 


This  is  an  alternate  form  of  the  Richards  equation  (Richards,  1931;  Jury  et  al.,  1991) 
using  effective  saturation  (Se)  as  the  dependent  variable. 

Assuming  that  air  exists  everywhere  at  atmospheric  pressure,  the  air  pressure 
head  is  zero  and  by  equation  (2-8): 


hc  =  ~K 


(2-21) 
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2.3.5  Zero-Viscosity  Model 


Equations  for  the  zero-viscosity  model  are  derived  by  evaluating  the  two- 
fluid  flow  equations  in  the  limit  as  the  air  viscosity  (  pa  )  goes  to  zero  (McWhorter  and 
Marinelli,  1996).  Evaluating  this  limit  for  equations  2-1 5b  through  2-1 5d,  we  obtain: 


fv(Se)  =  lim(Ma-0)  [f(Sa)]  = 


0 

1 


for  Sa-c1 
for  Sa= 1 


(2-22b) 


Dv(Se)  =  lim(Ma-0)  [D(Se)]  = 


kpwg  kJS9) 
Mw  C(Se) 

0 


for  Se<1 
for  Se  =  1 


=  lim(Ma-0)  [ E(Se) ]  = 


kPw9 

Mw 


U  se) 


for  Se*<  1 
for  Se  =  1 


(2-22c) 


(2-22d) 


In  each  case,  a  discontinuity  exists  between  values  associated  with  Se  =  1  and  values 
for  Se-1 .  Because  kfa(1)  =  0,  the  two-fluid  flow  equations  indicate  that  the  above 
values  given  for  Se  =  1  apply  for  any  air  viscosity.  Furthermore,  the  equations  for  Se  -c  1 
in  (2-22c)  and  (2-22d)  are  identical  to  those  defining  DR  and  ER  in  the  Richards  flux 
equation  (equation  2-19).  With  these  considerations,  the  flux  equation  for  the  zero- 
viscosity  model  is  written  as: 


q 


W 


q» 


for  Se  <  1 
for  Se  =  1 


(2-22a) 
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and  the  PDE’s  for  the  zero-viscosity  model  are: 


'dE^as. 
,  dSeJ  ax 


<j>(i-sr)^i  and 


ax 


0  for  Se<1 
for  Se=1 


(2-23) 


For  Se  <  1 ,  the  zero-viscosity  model  requires  that  two  differential  equations  be  solved 
simultaneously.  Note  that  equation  (2-23)  allows  for  the  possibility  of  a  saturation 
discontinuity  at  Se  =  1. 


2.3.6  Discussion 

Since  the  Richards  equation  is  sometimes  interpreted  to  imply  that  the  air 
viscosity  is  zero,  one  might  expect  the  Richards  and  zero-viscosity  models  to  be 
equivalent.  Inspection  of  the  flux  equations  and  PDE’s  for  both  models  (equations  2-19, 
2-20, 2-22,  and  2-23)  indicate  that  the  models  are  not  identical.  Thus,  the  traditional 
Richards  formulation  for  partially  saturated  flow  is  not,  strictly  speaking,  recovered  from 
the  two-fluid  flow  formulation  in  the  limit  of  zero  air  viscosity. 

Insight  into  the  differences  between  the  three  models  can  be  gained  by 
evaluating  equation  2-22b  for  the  function  f(Se)  at  different  water-air  viscosity  ratios 
(Mw/Ma).  0n  Figure  2-3,  the  function  is  plotted  for  the  typical  water-air  viscosity  ratio  of 
55.6  and  for  larger  values  (i.e.,  for  relatively  smaller  air  viscosities).  Note  that  f(1)  =  1 
for  any  viscosity  ratio.  This  feature  is  reasonable,  because  at  full  saturation,  there  is  no 
air  flow  since  kra(1)  =  0  and  water  flux  must  therefore  be  equal  to  the  total  fluid  flux.  In 
the  limit,  as  the  viscosity  ratio  goes  to  infinity  (air  viscosity  goes  to  zero),  f(1)  remains 
equal  to  unity,  but  a  discontinuity  develops  because  f  is  equal  to  zero  for  all  effective 
saturations  less  than  one.  This  is  the  zero-viscosity  model  definition  of  the  function  f  as 
described  by  equation  (2-22b).  In  the  Richards  model,  the  function  f  is  treated  as  equal 
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Figure  2-3.  Function  f(Se)  For  X  =  2. 


to  zero  for  all  effective  saturations,  which  is  why  the  function  does  not  appear  in  the 
Richards  equations. 

Inspection  of  flux  equation  (2-1 5a)  indicates  that  for  general  two-fluid  flow, 
knowledge  of  the  total  fluid  flux  (qt)  is  required  in  order  to  preserve  fluid  continuity.  As 
the  viscosity  ratio  goes  to  infinity,  equation  (2-22a)  for  the  zero-viscosity  model  shows 
that  knowledge  of  the  total  fluid  flux  is  needed  to  preserve  fluid  continuity  at  Se  =  1.  In 
the  Richards  equations,  knowledge  of  the  total  fluid  flux  is  lost,  because  the  flow 
dynamics  of  air  is  not  considered. 

Figures  2-4  and  2-5  are  plots  of  the  dimensionless  forms  of  the  diffusivity  [D(Se)] 
and  gravity  [E(Se)]  functions  for  different  values  of  the  water-air  viscosity  ratio.  It  is 
observed  that  D(1)  =  0  and  E(1)  =  0  for  all  viscosity  ratios.  These  values  are  retained  in 
the  limit  as  the  viscosity  ratio  goes  to  infinity,  which  represents  the  zero-viscosity  model. 
In  this  case,  each  function  exhibits  a  discontinuity  between  the  associated  values  for 
Se- 1  and  S0  =  1.  For  Se  ■<  1 ,  the  diffusivity  and  gravity  functions  for  the  Richards 
model  are  identical  to  those  of  the  zero-viscosity  model.  At  Se  =  1,  however,  the 
Richards  functions  are  continuous  and  do  not  exhibit  the  discontinuities  associated  with 
the  zero-viscosity  model. 

2.4  VERTICAL  DRAINAGE  IN  A  HOMOGENEOUS  SOIL 

To  further  investigate  the  failure  of  the  Richards  equation  and  the  overall  effect 
of  air  viscosity  on  water  drainage,  we  apply  the  two-fluid  and  Richards  models  to  the 
physical  problem  shown  on  Figure  2-1.  Due  to  the  discontinuities  in  functions 
associated  with  the  zero-viscosity  model,  it  cannot  be  applied  directly.  However,  the 
effect  of  reduced  air  viscosity  can  be  investigated  by  examining  the  behavior  of  the  two- 
fluid  model  as  the  air  viscosity  is  progressively  decreased  to  very  small  values. 
Simulations  are  performed  for  the  case  where  B  ■<  hd,  which  can  be  evaluated  using 
the  Richards  model,  and  for  B  >-  hd  which  is  not  solvable  by  the  Richards  model. 
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Figure  2-4.  Dimensionless  Diffusivity  Function  for  A  =  2  and  =  15  cm  (water). 
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Figure  2-5.  Dimensionless  Gravity  Function  for  A  =  2. 


The  initial-boundary-value  problem  described  on  Figure  2-1  can  be  simulated  for 
both  the  two-fluid  and  Richards  models  using  semi-analytical  solutions  described  in 
Chapter  1  of  this  report.  This  problem  is  described  as  vertical,  one-dimensional, 
unidirectional  displacement  of  a  wetting  fluid  (water)  by  a  nonwetting  fluid  (air)  in  a 
homogeneous  soil.  Note  that  if  (1)  the  fluids  and  porous  medium  are  assumed  to  be 
incompressible,  (2)  the  bottom  of  the  column  is  fully  saturated,  and  (3)  there  is  no  water 
flow  at  the  top  of  the  column,  the  following  relationship  holds  for  one-dimensional, 
unidirectional  displacement: 

q,(t)  =  qa(L,t)  =  qw(x*,t.)  =  qJO.t)  (2-24) 

From  a  physical  standpoint,  this  equation  implies  that  the  sides  of  the  column  are 
impermeable  and  that  air  access  into  the  system  only  takes  place  at  the  top  of  the 
column. 

First,  we  consider  drainage  in  a  100  cm  column,  composed  of  medium-grained 
sand,  where  water  pressure  head  at  the  base  is  maintained  at  zero  [h^O.t)  =  B  =  0]. 
The  soil  and  fluid  properties  used  in  this  simulation  are  summarized  in  Table  2-2. 

Based  on  solutions  to  the  two-fluid  model,  Figure  2-6a  shows  computed  saturation 
profiles  at  an  elapsed  time  of  98  seconds  for  the  typical  water-air  viscosity  ratio  of  55.6 
and  for  higher  values  (lower  air  viscosities).  For  the  typical  viscosity  ratio,  the  extent  of 
air  invasion  is  from  the  top  of  the  column  down  to  the  desaturation  point  located  at  x*  = 
92.2  cm.  As  the  viscosity  ratio  becomes  larger,  the  extent  of  air  invasion  increases. 

For  viscosity  ratios  of  10®  and  larger,  the  saturation  distribution  does  not  change  and  we 
therefore  consider  this  profile  to  approach  the  behavior  of  the  zero-viscosity  model.  As 
discussed  in  the  previous  section,  the  parameter  B  for  this  example  is  sufficiently  small 
that  evaluation  by  the  Richards  model  is  theoretically  possible.  As  shown  on  Figure  2- 
6a,  the  saturation  profile 
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Table  2-2.  Parameters  Used  in  Chapter  2. 


Pw 

1 

gem'3 

Mw 

0.01 

-1  -1 
g  cm  1  s  1 

Ma 

0.00018 

-1  -1 
g  cm  1  s  1 

(lower  values  used  in  some 

simulations) 

g 

980 

-2 

cm  s  * 

Medium  Sand 

Fine  Sand 

k 

10‘7 

2 

cm4 

10'8  cm2 

X 

3 

1.5 

hd 

15 

cm 

30  cm 

sr 

0.1 

0.2 

4> 

0.35 

0.35 
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a.  Saturation 


b.  Air  Pressure  Head 


EXPLANATION 
•  Desaturation  Point 

a  Mw/Mo  =  55-6 

b  =  103 

c  =  104 

d  =  ^106;  Richards  Model 

Figure  2-6.  Drainage  in  a  Homogeneous  Soil. 
hw(0,t)  >  -  hj  ;  t  =  98  seconds 
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determined  by  solution  of  the  Richards  equations  is  identical  to  that  predicted  by  the 
two-fluid  model  at  low  air  viscosity.  Note,  however,  that  the  shape  of  saturation  profile 
for  the  typical  viscosity  ratio  is  substantially  different  from  predictions  based  on  the 
Richards  model  and  the  two-fluid  model  with  low  air  viscosity. 

Figure  2-6b  shows  computed  air  pressure  head  distributions  associated  with  the 
saturation  profiles  on  Figure  2-6a.  For  the  typical  viscosity  ratio  of  55.6,  significant 
negative  pressure  heads  are  predicted  to  develop,  reaching  magnitudes  of  about  2.6  cm 
in  vicinity  of  the  desaturation  point.  As  expected,  when  the  viscosity  ratio  increases,  the 
air  pressure  heads  become  closer  to  zero  throughout  the  zone  of  air  invasion.  For 
viscosity  ratios  greater  than  106,  considered  for  this  example  to  approach  zero-viscosity 
behavior,  the  air  pressure  heads  are  essentially  equal  to  zero  throughout  the  zone  of  air 
invasion.  Since  the  Richards  model  assumes  zero  air  pressure,  Figure  2-6b  once  again 
indicates  that  the  predictions  of  the  two-fluid  model  at  low  air  viscosity  approach  those 
based  on  the  Richards  equation.  However,  results  of  the  two-fluid  analyses  indicate 
that  the  assumption  of  low  air  viscosity  is  not  necessarily  valid  for  real  air-water 
systems. 

Figures  2-7a  and  2-7b  show  saturation  and  air  pressure  head  distributions  at  an 
elapsed  time  of  1057  seconds  after  initiation  of  drainage.  For  the  typical  water-air 
viscosity  ratio  of  55.6,  air  invasion  extends  from  the  top  of  the  column  down  to  x*  =  47.2 
cm.  The  characteristics  of  these  plots  are  similar  to  those  shown  on  Figure  2-6,  with 
the  exception  that  (1)  the  extent  of  air  invasion  is  greater,  (2)  for  typical  air  viscosity,  the 
magnitude  of  air  pressure  head  is  greater  (more  negative)  at  the  desaturation  point,  and 
(3)  the  two-fluid  model  approaches  the  Richards  model  at  a  smaller  viscosity  ratio  of 
about  104.  For  this  example,  Figure  2-7  indicates  that  as  the  drainage  process 
proceeds,  the  effects  of  air  viscosity  remain  significant  both  on  the  shape  of  the 
saturation  profile  and  the  distribution  of  air  pressures  within  the  zone  of  air  invasion. 
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We  now  consider  simulation  of  drainage  in  the  same  soil  column  where  water 
pressure  head  at  the  base  [hw(0,t)]  is  maintained  at  -20  cm.  Since  this  value  is  smaller 
than  -hd,  the  magnitude  of  the  water  flux  will  exceed  the  saturated  hydraulic  conductivity 
(K)  and  the  Richards  model  cannot  be  used  to  simulate  this  physical  system.  Saturation 
profiles  at  an  elapsed  time  of  98  seconds  (Figure  2-8a  and  2-8b)  show  an  interesting 
feature  for  this  case  where  the  water  flux  is  greater  than  the  saturated  hydraulic 
conductivity  of  the  soil  [  |qw*  |  =  Jq^O.t)!  >-  K].  As  the  viscosity  ratio  increases,  the 
extent  of  air  invasion  becomes  greater  and  a  long  zone  with  very  high  water  saturations 
(low  air  saturations)  extends  progressively  farther  into  the  medium.  Note  that  the 
saturation  profile  for  the  viscosity  ratio  of  108  is  actually  plotted  at  a  time  of  64  seconds, 
because  at  98  seconds  it  is  predicted  that  the  zone  of  air  invasion  would  extend  all  the 
way  to  the  base  of  the  column.  Figure  2-8  illustrates  that  for  relatively  high  water 
drainage  rates,  desaturation  in  a  true  air-water  system  is  grossly  different  from  what 
would  be  predicted  by  a  system  with  very  low  air  viscosity. 

Consider  further,  the  associated  air  pressure  head  distributions  shown  on  Figure 
2-8c.  As  with  the  previous  example,  significantly  negative  air  pressure  heads  are 
predicted  for  the  typical  viscosity  ratio  of  55.6,  and  the  magnitudes  of  the  pressure 
heads  decrease  as  the  ratio  is  initially  increased  above  this  value.  However,  at  a  very 
large  viscosity  ratio  of  108  (very  small  air  viscosity),  the  magnitude  of  air  pressure  heads 
are  predicted  to  increase.  This  seemingly  counter-intuitive  result  can  be  explained  as 

follows.  Substituting  equation  (8)  into  (1 1)  and  taking  the  derivative  of  the  air  pressure 
head  with  respect  to  x  gives: 


dx 


nd  q-(1*m) 

A  e 


dS. 


dx 


(2-25) 
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As  discussed  previously,  at  very  large  viscosity  ratios,  the  lower  portion  of  the  soil 
column  develops  a  long,  partially  saturated  zone  where  Se  •*  1 ,  which  implies  that 
dSe/dx  -  0.  Inspection  of  equation  (2-25)  would  therefore  indicate  that  within  this 
zone: 


dh,  ^  dh^ 
dx  dx 


(2-26) 


Simulation  results  for  the  108  viscosity  ratio,  indicate  that  in  the  lower  portion  of  the 
saturation  profile,  that  is  from  the  desaturation  point  at  43.1  cm  up  to  50  cm,  the  air 
pressure  gradient  is  within  2  percent  of  the  water  pressure  gradient.  As  the  viscosity 
ratio  approaches  infinity,  we  expect  that  the  two  derivatives  would  become  identically 
equal. 

Finally,  for  the  physical  problem  shown  on  Figure  2-1,  we  consider  direct 
application  of  the  Richards  model  with  no  qualifications  regarding  air  access  into  the 
flow  system.  As  stated  previously,  a  basic  assumption  of  the  Richards  model  is  that  air 
is  freely  accessible  to  all  portions  of  the  flow  system,  regardless  of  the  saturation 
distribution.  For  the  lower  boundary  condition;  hw(0,t)  >■  -hd,  it  is  shown  on  Figures  2-6 
and  2-7  that  the  Richards  model  gives  identical  results  to  the  two-fluid  model  at  low  air 
viscosity.  For  this  boundary  condition,  if  the  access  of  low-viscosity  air  is  limited  to  the 
top  of  the  soil  column,  there  appears  to  be  essentially  free  air  access  to  all  portions  of 
the  column  undergoing  desaturation,  and  the  two  models  predict  the  same  drainage 
behavior.  However,  if  the  lower  boundary  condition  is  lTw(0,t)  <  -hd,  the  two  models 
are  not  equivalent  as  shown  diagrammatically  on  Figure  2-9.  Because  the  Richards 
model  assumes  that  the  air  pressure  head  is  zero,  this  boundary  condition  requires  that 
substantial  desaturation  occurs  immediately  at  the  base  of  the  column.  As  shown  on 
Figure  2-9a,  the  Richards  model  predicts  simultaneous  desaturation  at  both  the  top  and 
bottom  of  the  column,  with  the  central  portion  remaining  saturated.  Obviously,  this 
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Figure  2-9. 


a.  Richards  Model  b.  Two-Fluid 

Model 


Diagramatic  Comparison  of  the  Richards  Model  and  the  Two-Fluid 
Model  at  Low  Air  Viscosity. 

hw(0,t)  <  -  hd 


III  -  58 


situation  could  only  occur  if  air  inlets  are  provided  along  the  sides  of  the  column,  which 
would  violate  the  physical  problem  as  originally  stated.  For  low  air  viscosity,  Figure  2-9b 
shows  desaturation  along  the  entire  length  of  the  column,  but  with  saturation 
approaching  1  at  the  base.  This  is  equivalent  to  capillary  pressure  head  (hc)  being 
approximately  equal  to  the  displacement  pressure  head  (h^)  at  base  of  the  column, 
which  for  the  stated  water  pressure  boundary  condition,  requires  that  the  air  pressure 
head  is  less  than  zero.  As  shown  on  Figure  2-8,  under  these  conditions,  a  low  air 
viscosity  simulation  can  predict  air  pressure  variations  within  the  flow  system  even 
though  the  air  viscosity  is  very  small. 

2.5  VERTICAL  DRAINAGE  IN  A  LAYERED  SOIL 

At  face  value,  the  physical  problem  described  on  Figure  2-1  may  seem 
somewhat  artificial  because  of  the  negative  water  pressure  head  condition  imposed  at 
the  base  of  the  soil.  However,  it  can  be  shown  that  for  many  practical  drainage 
problems  in  layered  soils,  the  conditions  shown  on  Figure  2-1  can  occur  within  individual 
layers.  This  is  particularly  true  within  a  fine-grained  layer  that  overlies  a  coarser,  more 
permeable  layer.  Consider  the  layered-soil  drainage  problem  described  on  Figure  2-10 
for  the  case  of  real  air  viscosity.  The  properties  of  the  two  sands  that  comprise  the 
column  are  summarized  in  Table  2-2.  Initially,  the  column  is  fully  saturated  and 
hydrostatic,  with  the  water-table  (zero  water  pressure  head)  surface  located  at  the  top  of 
the  soil  column.  At  time  zero,  the  water  table  is  instantaneously  lowered  to  the  bottom 
of  the  column. 

Based  on  the  two-fluid  semi-analytical  solution  developed  in  Chapter  1 , 
computed  saturation  and  pressure  head  distributions  are  shown  on  Figure  2-1 1  at  an 
elapsed  time  of  3,360  seconds,  prior  to  arrival  of  the  desaturation  front  at  the  interface 
between  the  two  sands.  Of  particular  interest  is  that  significantly  negative  air  pressure 
heads  develop  within  the  zone  of  partial  saturation  and  that  no  desaturation  has 
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a.  Saturation 


b.  Pressure  Head 


h  (cm) 


EXPLANATION! 

•  Desaturation  Point 


Figure  2-11.  Drainage  in  a  Layered  Soil,  t  =  3360  seconds 


occurred  in  the  lower  layer.  Because  the  water  pressure  head  at  the  base  of  the  upper 
(fine  sand)  layer  is  less  than  the  negative  of  the  displacement  pressure  for  this  material, 
|qw*|  >  K^,  at  the  desaturation  point  and  the  Richards  formulation  cannot  be  applied  to 
this  drainage  problem.  This  occurs  even  though  the  water  pressure  head  at  the  base  of 
the  column  is  equal  to  zero.  Note  further  that  water  pressure  heads  in  the  underlying 
medium  sand  near  the  interface  are  well  below  the  value  at  which  desaturation  would 
occur  if  air  access  were  permitted  through  the  column  walls.  If  the  condition  of  no 
lateral  access  were  relaxed,  the  Richards  model  could  be  applied  and  it  would  predict 
that  the  lower  sand  would  have  long  since  desaturated.  However,  this  would  be  a 
departure  from  the  problem  statement  and  the  Richards  formulation  cannot  be  applied 
to  the  problem  as  posed. 

Upon  arrival  of  the  desaturation  front  at  the  material  interface,  the  lower  layer 
begins  to  drain.  Figure  2-12  shows  results  of  the  two-fluid  simulation  at  an  elapsed  time 
of  6,470  seconds.  Note  that  significantly  negative  air  pressures  persist  in  both  layers. 
The  existence  of  negative  air  pressure  heads  throughout  the  system  implies  that  for  a 
given  saturation  (given  capillary  pressure  head),  the  water  pressure  head  is  less  than 
that  which  would  occur  if  air  were  present  at  zero  pressure  head.  The  more  negative 
water  pressure  head  has  the  effect  of  decreasing  the  magnitude  of  the  drainage  rate, 
relative  to  the  rate  which  would  occur  assuming  zero  air  pressure. 

Predicted  drainage  rates  from  the  layered  soil  are  shown  on  Figure  2-13.  When 
the  desaturation  front  first  arrives  at  the  interface,  at  an  elapsed  time  of  about  4,600 
seconds,  there  occurs  a  transient  increase  in  the  magnitude  of  the  drainage  rate.  This 
is  caused  by  the  sudden  accessibility  of  air  to  the  lower  layer,  resulting  in  rapid 
desaturation  below  the  interface.  This  rapid  desaturation  is  manifested  as  an  increase 
in  the  magnitude  of  the  drainage  rate.  As  drainage  proceeds,  the  magnitude  once  again 
declines  and  eventually  approaches  zero  as  a  new  static  condition  is  approached.  This 
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a.  Saturation 


b.  Pressure  Head 


EXPLANATION 
•  Desaturation  Point 


Figure  2-12.  Drainage  in  a  Layered  Soil,  t  =  6470  seconds 
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Magnitude  of  Total  Flux:  qt  (cm/s) 


Note: 

Total  Flux  [qt(t)]  =  Air  (Inflow)  Flux  at  Top  of  Column  [qa(L,t)] 

=  Water  (Outflow)  Flux  at  Bottom  of  Column  [qw(0,t)] 


Figure  2-13.  Total  Volume  Flux  for  Drainage  in  a  Layered  Soil. 
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drainage  behavior  can  only  be  explained  by  considering  two-fluid  flow  with  finite  air 
viscosity. 

This  simulation  clearly  shows  that  for  drainage  in  a  layered  soil  column  with 
limited  access  of  viscous  air,  the  drainage  behavior  predicted  by  two-fluid  flow  theory  is 
grossly  different  than  what  would  be  predicted  if  one  assumes  either  zero  air  viscosity  or 
the  ability  of  air  to  freely  access  the  system  laterally.  For  the  case  of  viscous  air  with 
limited  access,  the  drainage  rates  are  lower  and  the  sequence  of  desaturation  from 
layer  to  layer  is  different.  The  Richards  model  provides  a  theoretical  solution  to  this 
problem  only  if  the  condition  of  no  lateral  air  access  is  relaxed. 

2.6  CONCLUSIONS 

From  the  standpoint  of  practical  application,  the  most  significant  conclusion  to 
be  drawn  from  this  study  is  that  the  hydraulics  of  air  flow  can  have  a  significant  impact 
on  the  nature  of  vertical  water  drainage.  The  effects  of  air  flow  appear  to  be  more 
pronounced  when  its  access  to  the  porous  medium  is  restricted  and  when  the  drainage 
process  occurs  in  layered  or  otherwise  heterogeneous  materials.  When  air  flow  is 
important,  rigorous  analysis  of  the  drainage  process  can  only  be  performed  by 
considering  two-fluid  (air-water)  flow,  with  air  possessing  a  finite  viscosity.  In  general, 
the  effects  of  viscous  air  and  limited  air  access  is  to  reduce  the  rate  of  drainage  and 
change  the  sequence  by  which  individual  layers  desaturate.  It  can  be  shown  that  for 
many  types  of  drainage  problems,  air  pressure  within  the  flow  system  is  reduced  to 
values  which  are  significantly  less  than  atmospheric. 

The  Richards  equation  has  traditionally  formed  the  basis  for  analysis  of  drainage 
in  air-water  systems.  When  air  access  is  limited,  the  Richards  model  can  lead  to 
inaccurate  simulations  of  drainage  behavior,  particularly  in  layered  soils.  Saturation 
distributions  and  drainage  rates  predicted  by  the  Richards  model  can  be  significantly 
different  from  the  real  behavior  of  field  soils  and  laboratory  test  columns.  For  a  certain 
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class  of  drainage  problems  with  limited  air  access,  the  Richards  model  does  not  allow  a 
theoretical  solution  to  be  obtained. 

Finally,  because  of  the  uniform  air  pressure  assumption  used  to  derive  the 
Richards  flux  equation,  the  Richards  model  could  be  interpreted  to  contain  an  implicit 
assumption  of  zero  air  viscosity.  In  our  opinion,  it  is  inaccurate  and  potentially 
misleading  to  interpret  the  Richards  formulation  as  describing  a  zero-air-viscosity  flow 
system.  It  is  more  accurate  to  view  Richards  equation  as  an  unlimited  air  access 
model,  where  air  is  supplied  to  the  system  at  zero  pressure  head.  There  are  many 
types  of  drainage  problems  where  the  Richards  model  will  provide  similar  results  as  a 
zero-air-viscosity  flow  model.  However,  for  certain  drainage  problems,  the  Richards 
model  does  not  predict  the  drainage  behavior  that  would  occur  in  a  true  zero-air- 
viscosity  system. 
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CHAPTER  3 


PERFORMANCE  AND  SIMULATION  OF  AN  LNAPL-WATER 
DRAINAGE  EXPERIMENT  IN  LAYERED  SOIL 


Abstract 

A  new  semi-analytical  solution  for  one-dimensional,  two-fluid  flow  is  used  to 
simulate  an  experiment  involving  the  drainage  of  LNAPL  and  water  in  a  laboratory 
column  containing  layered  sands.  Input  parameters  used  in  the  model  are  based  on 
independent  laboratory  measurements  of  the  fluid/soil  properties,  with  only  a  slight 
adjustment  in  the  hydraulic  conductivity  of  one  of  the  sands.  The  fit  between  model 
predictions  and  measured  drainage  rates  is  surprisingly  good,  considering  the  overall 
uncertainty  in  the  input  parameters  and  the  fact  that  the  model  is  not  extensively 
calibrated.  This  study  suggests  that  the  semi-analytical  solution  adequately  models 
two-fluid  drainage  in  layered  soils,  without  experiencing  the  types  of  stability  problems 
that  can  adversely  affect  the  performance  of  traditional  finite  difference  and  finite 
element  numerical  models. 

3.1  INTRODUCTION 

The  presence  of  Light  Nonaqueous  Phase  Liquids  (LNAPLs)  in  the  subsurface 
represents  a  source  of  ground-water  contamination  at  many  hazardous  waste  sites 
including  oil  refineries,  heavy  manufacturing  facilities,  petroleum  tank  farms,  and 
gasoline  service  stations.  Generally  distributed  near  the  water  table,  mobile  LNAPL  can 
migrate  vertically  when  the  water-table  elevation  changes.  Field  observations  have 
indicated  that  if  the  water  table  drops  below  its  previous  range  of  fluctuation,  LNAPL 
may  temporarily  or  permanently  disappear  from  monitoring  wells.  This  situation 
commonly  occurs  at  field  sites  after  initiation  of  a  ground  water  pump-and-treat  system, 
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which  can  lower  the  water  table  over  a  broad  area.  In  Marinelli  and  Dumford  (1996),  it 
is  proposed  that  the  above  observation  can  be  explained  by  the  process  of  slow 
downward  migration  of  LNAPL  above  the  (falling)  water  table.  It  is  further  considered 
that  the  LNAPL  migration  process  may  be  relatively  slower  and  exhibit  complex 
behaviors  when  it  takes  place  in  layered  or  otherwise  heterogeneous  soils. 

Simulation  of  vertical,  one-dimensional,  drainage  of  LNAPL  and  water  in  layered 
soils  has  been  difficult  using  traditional  finite  difference  and  finite  element  numerical 
models.  Due  to  the  highly  nonlinear  governing  equations,  these  modeling  approaches 
can  exhibit  numerical  instabilities  and  convergence  problems,  particularly  when  sharp 
saturation  fronts  move  across  material  interfaces.  Analytical  solutions  for  immiscible 
flow  developed  by  McWhorter  and  Sunada  (1990)  are  not  applicable  to  vertical  flow 
problems  and  do  not  explicitly  consider  layered  materials.  In  Chapter  1 ,  a  semi- 
analytical  solution  is  described  for  simulating  vertical,  one-dimensional,  two-fluid  flow  in 
layered  media.  The  two-fluid  semi-analytical  solution  has  been  successfully  applied  to 
layered  materials,  without  experiencing  the  types  of  stability  problems  that  can  adversely 
affect  the  performance  of  traditional  numerical  models. 

This  chapter  describes  a  simple  laboratory  drainage  experiment  where  water,  in 
the  absence  of  air,  is  displaced  by  LNAPL  in  a  layered  soil.  The  experiment  considers 
the  downward  migration  of  LNAPL,  below  a  shallow  ponding  condition,  due  to  a  rapid 
lowering  of  the  water  table.  The  experimental  setup  is  such  that  the  only  fluids  in  the 
system  are  LNAPL  (nonwetting  fluid)  and  water  (wetting  fluid).  The  observed  behavior 
of  this  drainage  experiment  is  compared  with  the  simulated  results  of  the  two-fluid 
semi-analytical  solution.  Independently  measured  properties  of  the  fluids  and  porous 
materials  are  used  as  input  to  the  semi-analytical  solution.  In  order  for  the 
comparison  to  be  objective,  the 
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mathematical  solution  is  not  extensively  calibrated  in  an  attempt  to  force  a  match 
between  the  simulation  results  and  laboratory  data. 

3.2  LABORATORY  EXPERIMENT 

The  laboratory  setup  used  for  the  LNAPL-water  drainage  experiment  is  shown 
on  Figure  3-1.  From  bottom  to  top,  the  laboratory  column  was  packed  with  40  cm  of 
100-mesh  silica  sand,  30  cm  of  40-60-mesh  silica  sand,  and  30  cm  of  100-mesh  silica 
sand,  all  manufactured  by  Colorado  Silica  Sands,  Inc.,  Colorado  Springs,  Colorado. 

The  sands  were  placed  by  gravity  raining  using  equipment  and  procedures  described  in 
Rad  and  Tumay  (1987).  After  sand  placement,  the  column  was  vacuum  saturated  with 
de-aired  water  and  left  overnight  so  that  any  remaining  air  would  dissolve  into  the  water 
phase.  The  ball  valve  was  then  opened  and  water  was  drained  from  the  column  until 
the  free  water  surface  was  just  above  the  perforated  plate.  A  Marriott  siphon  containing 
Soltrol  220  was  installed  at  the  top  of  the  column.  Soltrol  220  (Phillips  Petroleum 
Company,  Bartlesville,  Oklahoma)  is  a  very  stable  LNAPL  that  is  immiscible  with  water 
and  maintains  its  fluid  properties  when  in  contact  with  water.  The  Marriott  was 
positioned  so  that  a  Soltrol  pond  about  5  cm  deep  was  maintained  at  the  top  of  the 
column. 

At  the  beginning  of  the  experiment,  the  ball  valve  was  opened  allowing  water  to 
discharge  from  the  bottom  of  the  column  which  caused  the  LNAPL-water  interface  at 
the  top  of  the  column  to  approach  the  upper  surface  of  the  sand.  Zero  time  was 
specified  when  the  interface  reached  the  top  of  the  upper  sand  layer.  LNAPL  flowing 
into  the  top  of  the  sand  column  displaced  water  which  was  draining  downward  and 
discharging  from  the  base  of  the  column.  A  graduate  cylinder  was  used  to  measure  the 
volumes  of  water  discharging  from  the  base  of  the  column.  Water  drainage  was 
monitored  at  ten  minute  intervals  for  approximately  12  hours. 
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Marriott  Siphon 
Containing  Soltrol  220 


Figure  3-1.  Laboratory  Setup  for  the  LNAPL-Water  Drainage  Experiment. 


Water  discharge  volumes,  measured  at  periodic  time  intervals,  were  used  to 
computed  the  water  discharge  flow  rate  as  a  function  of  time.  Water  drainage  flux 
(darcy  velocity)  was  computed  as  this  flow  rate  divided  by  the  internal  cross-sectional 
area  of  the  column.  Assuming  that  the  fluids  and  porous  medium  are  incompressible, 
the  water  flux  at  the  base  of  the  column  should  be  approximately  equal  to  the  LNAPL 
flux  at  the  top  of  the  column.  Water  drainage  flux  over  time  is  shown  graphically  on 
Figure  3-2.  Of  interest  is  the  decreasing  water  flux  during  the  first  10,000  seconds,  the 
uniform  and  slowly  decreasing  flux  between  10,000  and  34,000  seconds,  and  the  sharp 
decrease  in  flux  observed  at  about  34,500  seconds. 

3.3  FLUID  AND  SOIL  PROPERTIES 

Properties  of  the  fluids  and  soils  used  in  this  experiment  have  been  measured  by 
Miller  and  Dumford  (1996)  and  are  summarized  in  Table  3-1  a, b.  Dry  bulk  densities  of 
the  sand  layers  in  the  experimental  column  (Table  3-1  c)  were  determined  by  weighing 
the  column  before  and  after  the  placement  of  each  layer.  The  average  porosity  of  the 
sand  layers  was  calculated  based  column  weights  before  and  after  saturation. 

3.4  MATHEMATICAL  SIMULATION  OF  THE  EXPERIMENT 

The  governing  flow  equations  for  one-dimensional,  incompressible,  two-fluid  flow 
are  discussed  in  detail  in  Chapter  1.  The  initial  condition  and  boundary  conditions 
describing  the  laboratory  experiment  are  shown  on  Figure  3-3.  Initially  the  layered 
column  is  completely  saturated  with  water.  The  boundary  conditions  at  the  top  of  the 
column  are  (1)  no-flow  of  water  and  (2)  a  nonwetting  fluid  (LNAPL)  pressure  head  is 
equal  to  3.98  cm  (water).  The  latter  condition  is  equivalent  to  a  Soltrol  pond  height  of 
4.8  cm  above  the  top  of  the  sand.  At  the  base  of  the  column,  the  boundary  condition  is 
a  water  pressure  head  equal  to  -2.3  cm,  which  is  the  elevation  difference  between  the 
discharge  pipe  and  the  base  of  the  fine  sand.  It  is  assumed  that  there  are  no  head 
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Water  Drainage  Flu 


Figure  3-2.  Drainage  Experiment  -  Measured  Drainage  Flux. 


Table  3-1.  Properties  of  Fluids  and  Soils  Used  in  the  Drainage  Experiment. 


a.  Fluids  (a) 

Density  (b) 

Viscosity  (b) 

Surface  Tension 
Interfacial  Tension 

b.  Soils  (a) 

Sand 

Description 

Mean  Grain  Diameter,  d50 

d10 

d90 

Uniformity  Coefficient  (dgQ/d^) 

Brooks-Corey  Air-Water 
Displacement  Pressure  Head 

Brooks-Corey  Pore  Size 
Distribution  Index  (b) 

Irreducible  Saturation  (b) 

Average  Packed  Porosity 

Saturated  Hydraulic 
Conductivity  (b,c) 

c.  Soils  (Experimental  Column) 

Sand 

Dry  Bulk  Density 

Porosity  (Average  for  Both  Sands) 


Units 

Water 

Soltrol  220 

(g  cm'3) 

1.00 

0.83 

(g  cm'1  s'1) 

0.0095 

0.0365 

m4 

(dyne  cm  ') 

70 

28 

(dyne  cm"1) 

N/A 

42 

Units 

Fine  Sand 

Medium 

100-Mesh 

40-60-Mesh 

(mm) 

0.16 

0.27 

(mm) 

0.08 

0.25 

(mm) 

0.30 

0.30 

2.18 

1.12 

(cm  water) 

41 

18 

1.7 

3.4 

0.08 

0.05 

0.37 

0.37 

(cm  s  ') 

0.0020 

0.010 

Units 

Fine  Sand 

Medium 

(g  cm'3) 

1.62 

1.67 

0.36 


Notes: 

(a)  Values  reported  in  Miller  and  Dumford  (1996) 

(b)  Input  values  for  the  simulations. 

(c)  Second  simulation  used  hydraulic  conductivity  of  0.0018  cm/s  for  fine  sand. 
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Figure  3-3.  Initial-Boundary-Value  Problem  for  Simulation  of  the  Drainage 
Experiment. 
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losses  in  the  coarse  sand  filter  below  the  fine  sand.  The  above  initial-boundary-value 
problem  is  solved  using  the  semi-analytical  technique  described  in  Chapter  1. 

Equations  relating  capillary  pressure,  fluid  saturation,  and  relative  permeabilities 
are  based  on  the  Brooks  and  Corey  (1964, 1966)  parametric  model.  For  this  physical 
problem,  the  wetting  fluid  (water)  is  always  under  drainage  conditions.  Therefore,  the 
capillary  pressure  -  saturation  relationships  are  assumed  to  be  nonhysteretic  and  to 
follow  main  drainage  paths.  Table  3-1  provides  air-water  displacement  pressure  heads 
for  the  sands  used  in  the  column.  For  this  simulation,  we  require  an  estimate  of  the 
LNAPL-water  displacement  pressures.  These  estimates  are  obtained  by  scaling  the  air- 
water  displacement  pressures  as  follows: 


lw 


haw 

H 


glW 

\^i 


(3-1) 


where  hj*  and  h™  are  the  displacement  pressure  heads  for  LNAPL-water  and  air- 
water  systems,  respectively,  olw  is  the  LNAPL-water  interfacial  tension,  and  o8*  is  the 
air-water  surface  tension.  For  the  Brooks-Corey  saturation  model,  this  procedure  is 
equivalent  to  the  scaling  approach  described  by  Parker  et  at.  (1987),  although  they 
apply  this  procedure  only  to  the  van  Genuchten  (1980)  parametric  model.  Based  on 
equation  (3-1),  the  LNAPL-water  displacement  pressure  heads  are  computed  to  be  10.8 
and  24.6  cm  (water)  for  the  medium  and  fine  sands,  respectively. 

Based  on  the  two-fluid  semi-analytical  solution,  the  laboratory  experiment  has 
been  simulated  using  the  above  LNAPL-water  displacement  pressures  and  other  input 
parameters  identified  in  Table  3-1 .  On  Figure  3-4,  simulated  water  drainage  rates  from 
two  modeling  runs  are  compared  with  actual  drainage  rates  measured  during  the 
laboratory  experiment.  The  first  simulation  is  based  on  the  initial  best-estimate 
hydraulic  conductivity  of  0.002  cm/s  for  the  fine  sand. 
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Water  Drainage  Flux  (cm/s) 


Explanation 

*  Measured 

Predicted  ;  K  (fine  sand)  =  0.0020  cm/s 
■  Predicted  ;  K  (fine  sand)  =  0.00184  cm/s 

Figure  3-4.  Comparison  of  Predicted  and  Measured  Drainage  Flux. 
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It  is  observed  that  the  predicted  fluxes  are  consistently  higher  than  the  measured  fluxes. 
The  second  modeling  run  is  based  on  a  reduced  hydraulic  conductivity  of  0.0018  cm/s 
for  the  fine  sand.  This  decrease  in  conductivity  of  approximately  8  percent  is 
considered  to  be  within  the  uncertainty  range  associated  with  the  laboratory 
permeameter  test  used  by  Miller  and  Dumford  (1996)  to  measure  the  conductivity  of  the 
fine  sand.  The  slightly  lower  conductivity  may  be  reasonable  because  the  packed 
porosity  of  the  experimental  column  (36  percent)  was  less  than  the  packed 
permeameter  porosity  (37  percent).  The  second  simulation  gives  lower  fluxes  and 
provides  a  better  fit  to  the  measured  data. 

Based  on  the  second  model  simulation,  Figure  3-5  shows  predicted  saturation 
profiles  at  different  times  during  the  drainage  process.  At  a  time  of  3,300  seconds 
(Figure  3-5a),  the  desaturation  front  is  in  the  uppermost  fine  sand  layer,  Figure  3-5b 
shows  the  profile  at  11,570  seconds,  just  after  the  desaturation  front  has  encountered 
the  first  material  interface.  Note  that  a  saturation  discontinuity  exists  across  the 
interface  and  that  the  saturation  approaches  one  at  the  base  of  the  upper  fine  sand 
layer.  At  20,570  seconds  (Figure  3-5c),  the  desaturation  front  has  moved  downward 
into  the  middle  medium  sand  layer.  The  saturation  discontinuity  remains  at  the  first 
interface.  On  Figure  3-5d,  at  32,970  seconds,  the  desaturation  front  has  just  moved 
past  the  second  interface.  Saturation  discontinuities  are  shown  to  exists  at  both 
interfaces  and  the  saturation  remains  close  to  unity  at  the  base  of  the  upper  layer.  In 
addition,  the  middle  layer  is  mostly  desaturated.  Finally,  at  44,970  seconds,  the 
desaturation  front  is  predicted  to  move  well  into  the  lower  fine  sand  layer.  The 
saturation  discontinuities  remain  at  both  interfaces.  However,  some  desaturation  is 
shown  to  occur  at  the  base  of  the  upper  layer. 
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3.5  DISCUSSION 


The  comparison  of  simulated  results  and  laboratory  data  on  Figure  3-5  shows 
that  the  semi-analytical  solution  seems  to  capture  the  qualitative  behavior  of  the 
laboratory  experiment.  The  general  shapes  of  the  curves  are  similar,  including  the 
sharp  reduction  in  water  flux  which  is  predicted  to  occur  when  the  desaturation  front 
migrates  into  the  lower  layer.  In  general,  the  semi-analytical  solution  predicts  higher 
fluxes  compared  to  those  measured  during  the  laboratory  experiment. 

The  second  simulation  provides  a  better  fit  to  the  laboratory  data.  These 
simulated  results  are  very  close  to  the  laboratory  measurements  during  the  first  14,000 
seconds.  After  14,000  seconds,  the  simulation  shows  fluxes  that  decrease  very  slowly, 
while  the  laboratory  experiment  indicates  more  rapidly  decreasing  fluxes.  The  abrupt 
decrease  in  flux  is  predicted  by  the  model  to  occur  at  about  31 ,000  seconds,  while  the 
observed  decrease  occurs  at  about  34,500  seconds.  The  higher  simulated  fluxes 
appear  to  have  the  effect  of  decreasing  the  time  for  the  middle  medium  sand  layer  to 
desaturate  compared  to  the  experimental  data.  As  a  result,  the  model  predicts  that  the 
saturation  front  reaches  the  second  interface  sooner  than  what  apparently  occurred 
during  the  experiment.  At  late  times,  the  simulated  drainage  fluxes  are  similar  to,  but 
slightly  higher  than,  the  laboratory  fluxes. 

Note  that  the  second  simulation  is  based  on  independently  measured  fluid  and 
soil  properties,  with  only  a  slight  adjustment  in  the  assumed  hydraulic  conductivity  of  the 
fine  sand.  As  a  consequence,  this  simulation  does  not  represent  an  extensive 
calibration  effort  to  force  a  match  between  simulation  results  and  measured  laboratory 
data.  Considering  the  overall  uncertainty  of  the  fluid/soil  properties,  the  predictions  of 
the  semi-analytical  solution  are  surprisingly  good. 
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3.6  CONCLUSIONS 


This  study  provides  evidence  that  the  semi-analytical  solution  for  two-fluid  flow 
adequately  models  the  important  physical  processes  controlling  LNAPL-water  drainage  in 
layered  soils.  The  fit  between  model  predictions  and  measured  laboratory  data  is  quite 
good,  considenng  that  the  input  values  are  uncertain  and  that  the  model  was  not 
extensively  calibrated.  The  semi-analytical  solution  achieves  this  level  of  performance 
without  experiencing  the  types  of  stability  problems  that  can  adversely  affect  the 
performance  of  traditional  finite  difference  and  finite  element  numerical  models.  This 
study  suggests  that  the  semi-analytical  solution  developed  in  Chapter  1  may  be  a 

promising  technique  for  modeling  the  one-dimensional  flow  of  immiscible  fluids  in  layered 
porous  media. 
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SECTION  IV 


SECTION  IV 


Chapter  1 

Stochastic  Aggregation  Model  (SAM)  for  NAPL-water 
displacement  in  porous  media 


Abstract 

Characterization  of  NAPL  distribution  in  the  subsurface  is  a  necessary  part  of 
managing  environmental  risks  posed  by  NAPLs.  Modeling  NAPL-water  displacement 
processes  at  the  pore-scale  can  lead  to  greater  insight  and  understanding  of  lateral 
spreading,  the  extent  of  vertical  penetration,  and  the  final  fluid  distributions  of  NAPL  in 
aquifers.  A  stochastic  aggregation  model  (SAM),  based  on  diffusion  limited  aggregation 
(DLA),  has  been  developed  that  can  predict  the  configuration  of  the  front  produced  when 
a  NAPL  displaces  water  at  the  pore-scale  in  a  homogeneous  porous  medium.  The  model 
takes  into  account  the  essential  properties  governing  front  stability:  NAPL-water 
viscosity  differences,  NAPL-water  density  differences,  intrinsic  permeability  of  the 
porous  media,  flow  rate,  and  the  inclination  angle  of  the  porous  media  from  the 
horizontal.  These  properties  are  the  independent  variables  for  the  transition  number, 
which  is  derived  from  two-phase  flow  equations  and  uses  a  macroscopic  non- 
dimensionalization.  The  transition  number  serves  as  the  model’s  input. 

Using  the  continuum  approach  to  determine  fluid-fluid  displacements  at  the  pore- 
scale  is  computationally  prohibitive  and  involves  many  assumptions.  Our  approach 
utilizes  the  stochastic  nature  of  the  DLA  algorithm,  while  including  fluid  properties, 
generalized  porous  media  properties  and  macroscopic  flow  rates.  Due  to  the  simplicity  of 
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the  algorithm,  this  modeling  technique  can  very  quickly  produce  realistic  NAPL-water 
configurations  and  allows  for  a  Monte  Carlo  approach  to  be  utilized. 

The  model  simulations  in  this  paper  show  front  configurations  for  common 
NAPLs  (1,2  dichloroethane,  1,2  dibromoethane,  gasoline,  and  mineral  oil)  which  range 
from  fingered  to  stable  and  can  be  compared  via  their  fractal  dimension.  The  simulations 
show  the  NAPL-water  displacement  front  configuration  as  it  is  influenced  by  the 
competition  between  capillary,  viscous,  and  gravitational  forces  at  the  pore-scale. 

1.1  INTRODUCTION 

Dense  Non-Aqueous  Phase  Liquids  (DNAPLs)  are  liquids  with  higher  specific 
gravities  than  that  of  water.  Their  use  has  been  widespread  since  World  War  II.  Most 
DNAPLs  are  solvents,  but  there  are  also  some  liquids  used  in  chemical  manufacturing, 
organic  synthesis,  and  mineral  separation  that  are  classified  as  DNAPLs.  Much  of  the 
contamination  caused  by  DNAPLs  has  been  the  result  of  leaking  underground  storage 
tanks  and  various  other  spills  caused  by  the  electronic,  instrument  manufacturing  and 
aerospace  industries  (Pankow,  et  al.  1 996).  While  their  designation  as  non-aqueous 
implies  immiscibility  with  water,  DNAPLs  can  cause  significant  contamination.  The 
absolute  solubilities  with  water  are  low  for  most  DNAPLs,  but  their  solubility  to 
maximum  contaminant  level  (MCL)  ratios  are  often  high. 

Light  Non-Aqueous  Phase  Liquids  (LNAPLs)  are  liquids  with  a  specific  gravity 
less  than  that  of  water.  Most  LNAPLs  are  associated  with  the  production,  refining,  and 
distribution  of  petroleum  products.  When  LNAPLs  are  released  at  the  surface  they  tend 
to  migrate  downward  through  the  unsaturated  zone  to  the  water  table.  If  large  gradients 
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are  present,  they  may  penetrate  the  water  table.  Petroleum  LNAPLs  are  relatively 
biodegradable,  but  their  constituents  such  as  benzene,  toluene,  ethyl  benzene,  and  xylene 
can  cause  widespread  contamination  (Bedient,  1994). 

Managing  the  environmental  risks  posed  by  NAPLs  and  improving  technologies 
used  to  clean  up  contaminated  areas  depends  largely  on  understanding  and  predicting  the 
migration  of  NAPLs.  Successful  NAPL  remediation  depends  on  determining  the 
presence  and  configuration  of  the  NAPL  phase.  In  the  range  of  possible  flow  regimes  and 
fluid  properties,  NAPL  displacement  of  water  in  porous  media  can  produce  fronts  which 
are  stable,  highly  fingered  or  anything  in  between.  The  movement  of  the  front  is 
determined  by  the  competition  between  viscous,  capillary,  and  gravity  forces  at  the  pore 
level.  By  understanding  the  movement  of  the  NAPL  at  the  pore-scale  we  can  gain  insight 
into  the  distribution  of  NAPLs  in  the  subsurface.  Specifically,  site  characterization  and 
the  level  of  detail  needed  for  field  sampling,  estimation  of  the  potential  groundwater 
contamination,  and  retention  capacities  can  be  better  determined. 

Modeling  the  NAPL  distribution  as  it  displaces  water  at  the  pore  level,  given  the 
characteristics  of  the  NAPL  and  release  conditions,  was  the  focus  of  this  research. 

Instead  of  using  the  continuum  approach,  a  stochastic  modeling  technique  was  chosen 
which  utilizes  basic  instability  criteria  for  two  fluids.  The  advantage  of  this  technique  is 
the  speed  with  which  the  algorithm  can  produce  simulations.  The  following  objectives 
were  determined  to  be  important  in  creating  a  model  that  could  be  useful  for  practical 
applications. 
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1 .  The  model  should  be  a  computationally  simple  pore-scale  model. 

2.  The  model  should  be  applicable  over  a  wide  range  of  fluid  properties  and 
flow  conditions. 

3.  The  model  should  have  macroscopic  input  parameters 

Our  stochastic  aggregation  model  (SAM)  is  a  variant  of  Witten  and  Sander’s 
(1983)  Diffusion  Limited  Aggregation  (DLA)  model.  DLA  has  been  used  to  model  such 
unstable  processes  as  dielectric  breakdown  (Niemeyer,  et  al.  1984;  Wiesmann  and  Zeller, 
1986),  diffusion-controlled  polymerization  (Kaufman  et  al.,  1986),  chemical  dissolution 
processes  (Daccord,  1986),  solute  leaching  (Flury  and  Fluhler,  1995),  viscous  fingering 
(Maloy,  et  al.,  1985;  Chen  and  Wilkinson,  1985),  fluid-fluid  displacement  in  Hele  Shaw 
cells  (Kadanoff,  1985;  Liang,  1986);  and  fluid-fluid  flow  in  porous  media  (Paterson, 
1984;  Lenormand,  1987;  Kiriakidis,  et  al.,  1991;  Fernandez,  et  al.,  1991).  The  input  to 
the  model  is  the  transition  number,  T,  which  is  a  linear  combination  of  the  macroscopic 
bond  and  capillary  numbers.  It  is  equivalent  to  the  instability  criteria  of  Saffman  and 
Taylor  (1958)  and  Chouke,  et  al.  (1959).  The  transition  number  accounts  for  disparities 
in  the  viscosities  and  densities  of  the  two  fluids,  and  its  use  as  the  input  to  the  modified 
DLA  algorithm,  SAM,  is  unique. 

1.2  BACKGROUND 

1.2.1.  Stochastic  Modeling  Techniques 

Lenormand,  et  al.  (1988)  performed  micromodel  experiments  and  showed  that 
dramatically  different  front  configurations  occur  during  fluid-fluid  displacement  in 
porous  media  that  can  be  modeled  by  different  stochastic  techniques.  The  domains  of  the 
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stochastic  modeling  techniques  are  shown  in  the  phase  diagram  in  Figure  1  (Lenormand, 
1987).  The  phase  diagram  is  applicable  to  displacement  of  the  wetting  fluid  by  the  non¬ 
wetting  fluid  without  gravity  forces.  The  three  domains  labeled  as  DLA,  IP,  and  Plug 
Flow  correspond  to  displacement  of  a  viscous  fluid  by  a  nonviscous  fluid,  capillary 
displacement,  and  displacement  of  a  nonviscous  fluid  by  a  viscous  fluid,  respectively. 
Plug  flow  is  modeled  by  a  technique  called  anti-DLA.  The  areas  between  these  pure 
domains  are  governed  by  a  combination  of  the  applicable  forces.  The  limits  of  these 
domains  are  determined  by  an  evaluation  of  the  capillary  pressure  at  the  interface 
between  the  two  fluids  and  use  only  the  mobility  ratio  (M)  and  the  microscopic  capillary 
number  (Ca)  as  parameters.  The  mobility  ratio,  M,  is  defined  as  the  ratio  of  the  dynamic 
viscosity  of  the  displacing  fluid  to  the  dynamic  viscosity  of  the  displaced  fluid: 


Mx 

The  microscopic  capillary  number,  Ca  is  defined  as  the  ratio  between  the  viscous  forces 
and  the  capillary  forces  (Lake,  1989): 


Ca=- 


<Uj 

CTCOS0 


(1-2) 


where  q  is  the  microscopic  volumetric  flux,  //  is  the  dynamic  viscosity,  a  is  the 
interfacial  tension  between  the  two  fluids,  and  6  is  the  contact  angle. 

Continuum  models  are  considered  by  Lenormand  (1987),  to  be  in  the  domain 
labeled  stable  flow.  When  immiscible  fluid-fluid  displacement  is  not  stable,  the  phase 
diagram  offers  modeling  techniques  for  the  pure  regimes  of  capillary  and  viscous 
fingering  and  the  limits  for  these  pure  regimes.  To  model  fluid-fluid  displacements 
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Figure  1.  Lenormand’s  (1987)  phase  diagram  for  immiscible  fluid-fluid  displacement 
in  the  absence  of  gravity.  The  stochastic  modeling  techniques  corresponding  to  the  regimes 
for  which  their  use  is  suggested  is  shown  with  varying  microscopic  capillary  number  and 

mobility  ratio. 


between  these  regimes,  the  combined  forces  acting  to  stabilize  or  destabilize  the  front 
must  be  considered.  The  effects  of  gravity  are  also  important  (Wilkinson,  1984; 
Birovljev,  et  al.,  1991)  and  should  be  included  in  determining  the  NAPL-water  front 
configuration. 

We  will  consider  the  modeling  techniques  of  IP  and  DLA  before  discussing  how 
to  combine  the  viscous  and  capillary  effects  as  well  as  incorporate  gravity  forces  into  the 
current  model,  SAM. 

Percolation  modeling  (Broadbent  and  Hammersley,  1957)  takes  place  on  a  lattice 
of  specified  geometry  and  as  such  is  designated  as  a  network  model.  Fatt  (1956)  was  the 
first  to  use  network  models  for  investigation  of  porous  media.  The  application  of 
percolation  modeling  to  immiscible  fluid-fluid  displacement  requires  two  assumptions. 
The  first  is  that  the  interface  moves  as  a  series  of  discrete  jumps  by  individual  menisci. 
Surface  tension  effects  purely  determine  these  jumps.  The  second  assumption  maintains 
that  the  porous  media  is  connected,  and  in  some  way  random  (Wilkinson,  1984).  In 
standard  percolation  modeling,  the  probability  that  a  lattice  site  will  be  filled  is  equal  to 
the  ratio  of  the  number  of  sites  currently  occupied  to  the  total  number  of  sites.  For 
immiscible  fluid-fluid  displacement,  each  lattice  site  has  a  value  assigned  to  it  that  is 
representative  of  the  pore  size.  When  the  algorithm  begins  (for  drainage  of  the  wetting 
fluid)  all  of  the  sites  with  a  value  greater  than  some  r*  will  drain;  where  r*  is  related  to 
the  capillary  pressure  in  the  system  (Wilkinson,  1984).  This  technique  is  limited  by  its 
implication  that  the  occupation  of  one  lattice  site  is  independent  of  the  occupation  of  its 
neighbor.  The  result  is  a  lattice  with  many  unconnected  clusters  where  the  invading  fluid 
has  displaced  the  defending  fluid.  In  invasion  percolation  (IP),  the  method  is  the  similar. 
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The  difference  is  that  for  a  given  capillary  pressure,  only  the  lattice  sites  on  the  interface 
between  the  two  fluids  may  be  invaded  and  regions  where  the  wetting  fluid  is  surrounded 
by  the  non-wetting  fluid  remain  trapped  (Wilkinson  and  Willemsen,  1983).  This  results 
in  a  continuous  cluster  instead  of  many  disconnected  ones.  IP  models  have  been  utilized 
to  investigate  hysteresis,  boundary  effects,  ganglia  and  blob  structure,  and  the  effects  of 
contact  angle  on  capillary  pressure  saturation  curves  (Chatzis  and  Dullian,  1985;  Li,  et 
al.,  1986;  Wardlaw,  et  al.,  1987;  Jerauld  and  Salter,  1990;  Tsakiroglou  and  Payatakes, 
1990;  Ferrand  and  Celia,  1992;  Soil  and  Celia,  1993;  Lowry  and  Miller,  1995). 

Berko  witz  and  Balberg  (1993)  present  a  good  overview  of  percolation  modeling  as  it  is 
applied  to  groundwater  hydrology.  IP  models  are  only  valid  for  flow  dominated  by 
capillary  forces,  i.e.  where  the  Young-Laplace  equation  governs  the  movement  of  the 
interface  (Dullien,  1992).  An  example  of  an  aggregate  produced  by  IP  is  shown  in  Figure 
2. 


Figure  2.  An  aggregate  produced  by  invasion 
percolation  with  trapping  (Meakin,  1991) 
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1.2.2  DLA  Modeling 

Aggregate  growth  by  irreversible  particle  collection  is  a  common  phenomenon  in 
nature  (Sander,  1984).  The  displacement  of  one  fluid  by  another  in  porous  media  can  be 
represented  by  particle  aggregation  where  the  growing  aggregate  represents  the  invading 
fluid  and  the  space  surrounding  the  aggregate  represents  the  defending  fluid.  Several 
modeling  techniques  can  be  used  for  particle  aggregation.  These  include  IP  modeling 
(Wilkinson  and  Willemsen,  1983),  Eden  modeling  (Sander  1984),  ballistic  aggregation 
(Sander,  1984),  and  DLA  (Witten  and  Sander,  1983).  IP  modeling  was  discussed  in  the 
previous  section,  and  has  can  be  applied  to  immiscible  fluid  displacements  dominated  by 
capillary  forces.  To  be  able  to  model  fluid-fluid  displacements  where  any  combination  of 
viscous,  capillary,  or  gravity  forces  can  govern  the  stability  of  the  interface,  we  must  take 
a  look  at  the  other  aggregate  modeling  options.  Neither  Eden  modeling  nor  ballistic 
aggregation  can  describe  the  growth  instabilities  which  are  characteristic  of  many 
processes  including  fluid-fluid  displacement  with  unfavorable  mobility  ratios  (Sander, 
1986).  DLA  models,  however,  produce  aggregates  resulting  from  unstable  phenomenon 
using  a  stochastic  process.  Following  Witten  and  Sander  (1983),  numerous  other 
researchers  have  expanded  on  both  the  modeling  technique  and  the  applications  for  the 
model  (Meakin  and  Deutch,  1986;  Meakin,  1986;  and  Jullien,  et  al.,  1984).  Meakin 
(1988)  provides  a  thorough  review  of  DLA  as  a  model  for  kinetic  growth  phenomena. 

Diffusion-limited  aggregation  is  an  algorithm  whereby  an  initial  “seed”  particle  is 
placed  in  the  center  of  a  grid.  One  at  a  time,  particles  are  released  from  a  random 
position  on  the  boundary.  The  particles  move  from  location  to  location  on  the  grid  by 
taking  one-unit  steps  in  a  randomly  chosen  direction.  When  the  particle  reaches  a 
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location  adjacent  to  a  particle  which  is  part  of  the  aggregate,  the  walking  particle  becomes 
part  of  the  aggregate,  and  another  particle  is  released. .  A  DLA  aggregate  is  shown  in 
Figure  3.  Qualitatively,  the  aggregates  that  are  formed  in  this  manner  appear  disordered 
with  no  length  scale.  This  is  a  result  of  the  “screening”  effect  of  the  branches  (Sander, 
1986).  As  the  aggregate  grows,  it  becomes  more  and  more  difficult  for  a  particle  to 
traverse  the  narrow  areas  between  the  branches  without  coming  in  contact  with  and 
thereby  becoming  part  of  the  aggregate.  Consequently,  more  particles  will  come  in 
contact  with  the  aggregate  at  the  ends  of  the  branches. 


Figure  3.  An  aggregate  produced  by  Witten  and  Sander’s  (1983) 
original  diffusion  limited  aggregation  algorithm. 

The  probability  distribution  of  the  random  walkers,  u(r,t),  obeys  the  following 
equation  in  the  continuum  limit  (Sander,  1986): 


where  r  is  a  distance  from  the  aggregate,  t  is  time,  and  D  is  a  diffusion  constant.  The 


boundary  conditions  are:  u  equals  a  constant  as  r  goes  to  infinity,  and  w  =  0  on  the 
aggregate  surface.  For  a  steady  flux  of  walkers  from  a  source  that  is  far  away  from  the 
aggregate  surface,  the  equation  can  be  reduced  to  the  Laplace  equation  (Sander,  1986): 


V2u  =  0 


(1-4) 


1.2.3  Variations  on  DLA 

There  have  been  many  publications  that  have  explored  the  effects  of  specific 
variations  on  the  original  DLA  algorithm.  The  authors  of  this  paper  have  employed  some 
of  these  in  the  exploration  of  the  adaptation  of  DLA  to  fluid-fluid  displacement  in  porous 
media  with  varying  flow  rates  and  fluid  properties. 

To  shorten  the  time  it  takes  walkers  to  reach  the  aggregate  from  very  far  away, 
Meakin  (1988)  allowed  particles  to  take  larger  steps  if  they  were  at  a  large  distance  from 
the  occupied  sites.  This  improved  the  efficiency  of  the  algorithm  and  allowed  the 
aggregates  to  be  constructed  faster  while  maintaining  the  same  general  morphology. 

Jullien,  et  al.  (1984)  were  the  first  to  use  a  seed  line  instead  of  a  seed  particle.  In 
lieu  of  sending  walkers  from  a  random  point  on  a  boundary  circle,  they  sent  walkers  from 
a  randomly  chosen  location  on  a  finite  width,  semi-infinite  “ribbon”  towards  a  finite  line 
of  seed  particles.  This  adaptation  produces  aggregate  formations  that  are  similar  to  flow 
paths  made  by  a  non-wetting  fluid  displacing  water  in  a  Hele-Shaw  cell.  A  DLA 
simulation  using  a  line  seed  is  shown  in  Figure  4. 

The  third  variation  that  we  will  present  is  the  sticking  probability.  When  a  two- 
dimensional  circular  object  is  related  to  its  mass  by  the  second  power  of  the  radius  we 
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call  that  object  compact.  If  a  two-dimensional  circular  object  such  as  a  DLA  aggregate  is 
irregular,  we  can  take  the  radius  of  gyration  to  determine  the  mass  density  for  that  object 
in  two  dimensions.  If  a  number  between  one  and  two  relates  the  radius  of  gyration  and 
the  mass,  we  say  that  the  object  is  fractal  (Ferer,  et  al.,  1995).  The  number  relating  the 
radius  of  gyration  and  the  mass  is  termed  the  fractal  dimension,  D,  and  can  be  computed 
by  several  different  techniques. 


This  notion  of  the  area  covered  by  an  object  in  two  dimensions  can  be  thought  of 
as  a  two  dimensional  density  of  the  object.  We  can  alter  the  density  of  the  aggregates 
created  in  a  DLA  simulation  by  changing  the  sticking  probability.  The  sticking 
probability  is  the  probability  that  a  random  walker  that  steps  onto  a  location  adjacent  to 
an  occupied  location  will  “stick”  and  become  part  of  the  aggregate.  The  sticking 
probability  essentially  determines  how  susceptible  the  interface  of  the  aggregate  is  to 
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rearrangements.  When  the  sticking  probability  is  unity,  we  observe  classic  DLA 
formations.  As  the  sticking  probability  is  reduced,  the  aggregate  becomes  more  compact 
which  results  in  a  thickening  of  the  aggregate’s  branches.  In  the  limit  as  the  sticking 
probability  goes  to  zero,  the  aggregate  will  qualitatively  approach  a  stable  front  in  the 
analogy  of  fluid-fluid  displacement  in  porous  media.  It  was  found  in  the  current  study 
that  this  transition  from  fractal  to  compact  configurations  follows  a  power  law 
relationship. 

1.2.4  Applications  of  DLA 

The  use  of  DLA  to  model  unstable  processes  which  have  similar  universal 
properties  has  been  investigated  (Sander,  1986).  These  include  diffusion-controlled 
polymerization  (Kaufman  et  al.,  1986),  chemical  dissolution  processes  (Daccord,  1986), 
dielectric  breakdown  (Niemeyer  et  al.,  1984;  Meakin,  1986),  solute  leaching  (Flury  and 
Fluhler,  1995),  viscous  fingering  (Maloy  et  al.,  1985;  Chen  and  Wilkinson,  1985),  and 
fluid-fluid  displacement  in  Hele  Shaw  cells  (Kadanoff,  1985;  Liang,  1986).  In  its  first 
application  to  fluid-fluid  flow  in  porous  media,  Paterson  (1984)  proposed  a  mathematical 
analogy  between  flow  in  porous  media  and  DLA.  The  macroscopic  flow  through  porous 
media  is  governed  by  Darcy’s  law: 

q  -  —  (Vp  -  pgS/z)  (1-5) 

For  an  incompressible  fluid,  the  divergence  of  the  flux,  q,  is  zero  and  in  the  absence  of 

I 

gravity  we  have  the  Laplace  equation  for  the  pressure  in  the  more  viscous  fluid: 
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(1-6) 


V2jp  =  0 

According  to  Paterson  (1984),  this  equation  is  analogous  to  Equation  (1-14)  when  the 
pressure  drop  in  the  invading,  less  viscous  fluid  is  zero.  Chen  and  Wilkinson  (1985) 
agreed  with  this  analogy  and  hypothesized  that  the  randomness  of  a  statistically 
homogeneous  porous  media  on  the  pore-scale  was  modeled  by  the  stochastic  nature  of  the 
random  walk.  The  authors  who  have  previously  performed  experiments  to  determine  if 
the  Laplace  analogy  between  DLA  and  flow  through  porous  media  is  correct  use 
Equation  (1-6)  as  the  basis  of  the  relationship  between  the  DLA  model  and  the  fluid-fluid 
displacement  in  porous  media  (Paterson,  1984;  Chen  and  Wilkinson,  1985;  Maloy,  et  al., 
1985;  Lenormand,  1987;  Fernandez,  et  al.,  1991;  Kiriakidis,  et  al.,  1991). 

The  history  of  using  DLA  to  model  fluid-fluid  displacement  began  with 
comparing  Hele-Shaw  experiments  to  DLA  simulations.  By  applying  the  fluid 
incompressibility  condition  to  the  flow  equation  for  a  viscous  fluid  being  displaced  in  a 
Hele-Shaw  cell,  Equation  (1-5)  is  obtained.  Kadanoff  (1985)  and  Liang  (1986)  were  the 
first  to  use  DLA  to  model  air  displacing  water  in  a  Hele-Shaw  cell.  Both  Kadanoff 
(1985)  and  Liang  (1986)  modified  the  original  DLA  algorithm  to  include  surface  tension 
effects.  This  essentially  made  the  fingers  look  more  rounded  and  less  dendritic  at  the 
ends.  Nittmann,  et  al.  (1985)  and  Daccord,  et  al.  (1986)  carried  out  radial  Hele-Shaw 
experiments  with  water  displacing  a  high  viscosity,  non-Newtonian  polymer  and 
compared  the  resulting  fluid  configurations  to  simulations  from  an  original  DLA  model. 
Their  simulations  very  closely  match  their  experimental  results  and  they  concluded  that 
DLA  quantitatively  modeled  the  experimental  fluid-fluid  configuration  as  measured  by 
the  fractal  dimension  of  each.  It  was  later  postulated  by  Maloy,  et  al.  (1985)  that  the 
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DLA  algorithm  without  variations  had  successfully  modeled  the  radial  Hele-Shaw  fluid- 
fluid  displacements  of  Nittmann,  et  al.  (1985)  because  the  required  randomness  was  due 
to  viscosity  fluctuations  in  the  more  viscous,  non-Newtonian  fluid.  Sherwood  (1987) 
used  the  mathematical  analogy  to  DLA  (the  Laplace  equation  for  pressure)  to  obtain  front 
configurations  for  the  motion  of  a  fluid  through  porous  media.  The  examples  generated 
by  his  model  are  remarkably  similar  to  the  examples  generated  by  a  DLA  algorithm. 

Maloy,  et  al.  (1985),  Lenormand  (1988),  Fernandez,  et  al.  (1991),  and  Kiriakidis, 
et  al.  (1991)  have  all  explored  the  use  of  DLA  modeling  for  fluid-fluid  displacement  in 
porous  media  using  Paterson’s  original  analogy. 

Maloy,  et  al.  (1985)  conducted  experiments  whereby  epoxy  was  displaced  by  air 
in  a  radial,  two-dimensional  porous  medium.  The  resulting  finger  structures  are  tree-like 
and  have  the  same  fractal  dimension  as  aggregates  produced  by  DLA.  In  their 
determination,  when  the  flow  rate  is  high  (more  dominant  viscous  forces),  the  interface 
advances  on  the  pore  level  by  selecting  the  most  favorable  pore.  They  (Maloy,  et  al., 
1985)  conclude  that  DLA  more  appropriately  models  this  behavior  because  the  aggregate 
growth  is  controlled  by  both  the  Laplace  equation  and  the  discrete  random  advance  on  the 
pore  level.  In  IP  modeling,  the  movement  is  only  governed  by  the  pore  or  pore  throat 
size,  and  the  pressure  gradient  is  neglected. 

Lenormand  (1988),  as  presented  earlier,  suggested  the  use  of  different  stochastic 
models  with  differing  front  configurations,  and  gave  limits  for  these  models  in  terms  of 
the  microscopic  capillary  number  and  mobility  ratio.  Using  the  limits  calculated  by 
Lenormand  (1988),  Kiriakidis,  et  al.  (1991)  formulated  an  equation  for  the  sticking 
probability  for  yet  another  variation  on  DLA  consisting  of  two  types  of  random  walkers: 


wetting  and  non-wetting.  Gravity  was  not  considered.  The  model  simulations  generated 
by  Kiriakidis,  et  al.  (1991)  showed  that  DLA  could  be  used  to  produce  front 
configurations  ranging  from  fingered  to  stable.  These  simulations  appear  to  be 
qualitatively  representative  of  fluid-fluid  displacement  in  porous  media.  However,  the 
fluid  properties  associated  with  the  simulations  are  not  representative  of  common  fluids. 
The  model  simulations  were  generated  for  water  being  displaced  by  non-wetting  fluids 
with  dynamic  viscosities  of  7.6e-5  cp  and  13  cp.  With  the  exception  of  PCB’s,  the 
majority  of  dynamic  viscosities  for  NAPLs  lie  in  the  range  of  0.3  cp  to  8  cp  (Mercer  and 
Cohen,  1990). 

Fernandez,  et  al.  (1991)  used  a  DLA  model  with  a  releasing  probability  (converse 
of  a  sticking  probability)  to  show  that  there  is  a  crossover  between  DLA  and  IP  as  the 
releasing  probability  increases.  This  is  equivalent  to  saying  the  NAPL-water  front,  as 
modeled  by  DLA,  becomes  more  stable  as  the  sticking  probability  decreases.  The 
correlation  between  the  releasing  probability  and  the  physics  of  the  fluid-fluid 
displacement  does  not  include  gravity  and  is  related  to  microscopic  properties. 

The  previous  work  devoted  to  the  application  of  DLA  to  fluid-fluid  displacement 
in  porous  media  has  shown  that  DLA  simulations  can  be  produced  which  are  both 
qualitatively  and  quantitatively  (via  comparison  of  the  fractal  dimension)  similar  to 
experimental  fluid-fluid  displacements.  In  addition,  the  sticking  probability  was 
identified  by  both  Kiriakidis  (1991)  and  Fernandez  (1991)  to  be  the  parameter  by  which 
the  movement  of  the  interface  can  be  related  to  the  DLA  model. 
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1.2.5  Macroscopic  Dimensionless  Numbers 

Macroscopic  dimensionless  numbers  were  utilized  as  model  input  for  SAM. 
Hilfer  and  0ren  (1996)  re-examined  the  multiphase  flow  equations  and  looked  at  the 
definition  of  non-dimensional  variables  for  the  microscale  and  macroscale.  When  using 
dimensional  analysis  for  the  Navier  Stokes  equation,  the  normalized  pressure  field  is 
traditionally  defined  as: 


p  =  MLp  o-?) 

Jc 

where  the  hat  indicates  the  normalized  pressure.  This  definition  neglects  the  potential 
dependence  on  saturation  history,  and  bases  macroscopic  pressures  on  macroscopic 
viscous  pressure  effects.  Hilfer  and  0ren  have  instead  normalized  the  macroscopic 
capillary  pressure  by  the  displacement  pressure  from  the  capillary  pressure-saturation 
curve.  The  displacement  pressure  is  the  pressure  that  corresponds  to  the  inflection  point 
on  the  characteristic  curve.  By  using  this  normalization,  the  saturation  history  may  be 
accounted  for  as  the  displacement  pressure  on  the  drainage  or  imbibition  curve.  The 

resulting  macroscopic  capillary  number,  Ca ,  is  defined  as: 

Ca  =  ^-  O'8) 

tyd 

where  q  is  a  macroscopic  flux  for  the  phase  of  interest,  1  is  a  macroscopic  length,  and  k  is 
the  intrinsic  permeability.  Using  this  analysis,  the  macroscopic  gravity  number  and 
gravillary  number  may  also  be  defined.  The  gravity  number  is  the  ratio  of  viscous  forces 
to  gravity  forces  and  is  defined  as: 
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P  gk 

The  gravillary  number  is  the  ratio  of  gravity  forces  to  capillary  forces.  The  gravillary 
number  is  the  familiar  bond  number  (Lake,  1989)  when  the  difference  in  densities  is  used. 
The  gravillary  number  is  defined  as: 


d-10) 


1.3  THE  STOCHASTIC  AGGREGATION  MODEL  (SAM)  AND  THE 
TRANSITION  NUMBER 

SAM  is  a  modified  DLA  model  that  makes  use  of  some  of  the  variations 
mentioned  in  the  background  section  to  model  front  configurations  in  fluid-fluid 
displacements.  Specifically,  our  algorithm  models  the  displacement  of  water  by 
DNAPLs  which  are  both  more  and  less  viscous  than  water,  and  LNAPLs  which  are  both 
more  and  less  viscous  than  water.  In  addition,  we  look  at  the  flow  of  these  fluids  under 
the  conditions  of  increasing  flow  rates  and  in  a  gravity  field.  The  front  configurations 
produced  under  these  conditions  range  from  viscous  fingers,  to  percolation  clusters,  to 

stable  fronts.  The  parameters  of  the  model  are  macroscopic  dimensionless  numbers  and 
fluid  properties. 

Our  model  uses  Jullien  s  (1984)  approach  of  using  a  seed  line  where  the  particle 
walks  from  a  point  on  the  grid  far  from  the  seed,  in  a  region  of  finite  width.  A  seed  line 
in  our  case  represents  a  line  source  of  NAPL.  We  can  also  use  a  single  seed  particle  in 
the  middle  of  the  top  row  of  the  grid  to  emulate  a  NAPL  point  source.  To  decrease  the 
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time  it  takes  for  a  particle  to  walk  to  the  vicinity  of  the  aggregate,  we  used  Meakin’s 
(1988)  approach  of  allowing  the  walker  to  take  larger  steps  when  further  from  the 
aggregate.  If  the  particle  wanders  more  than  seven  rows  past  the  furthest  occupied  space, 
the  particle  is  removed  and  another  one  is  released.  The  final  variation  incorporated  into 
SAM  is  the  sticking  probability.  The  sticking  probability  is  the  only  input  into  the  model 
and  its  relationship  to  the  transition  number  will  be  discussed  in  the  following  section. 

For  the  simulations  presented  here,  a  square  network  of  10,000  elements  was 
used.  Each  element  represents  a  single  pore.  The  pore  size  may  be  approximated  via  the 
Kozeny-Carmen  equation  (Dullien,  1992)  given  the  value  for  intrinsic  permeability  used 
in  the  transition  number.  For  a  network  of  this  size,  one  simulation  runs  in  0.5  to  4 
minutes  on  an  IBM  RS6000  depending  on  the  sticking  probability. 

When  the  simulation  starts,  particles  are  released  from  a  random  location  opposite 
the  seed  point  or  line.  The  particle  then  has  an  equal  probability  of  moving  up,  down, 
right,  or  left.  If  the  particle  encounters  the  sides  of  the  network,  it  “bounces  off’  and 
continues  walking.  When  the  particle  reaches  a  location  adjacent  to  an  occupied  location, 
it  is  incorporated  into  the  aggregate  with  a  probability  equal  to  the  sticking  probability, 
and  a  new  particle  is  released.  The  simulations  presented  here  were  stopped  when  the 
bottom  of  the  grid  was  reached,  instead  of  prescribing  a  specific  number  of  particles  to 
form  an  aggregate.  Utilizing  continuity  in  the  model  is  currently  being  studied.  Standard 
IMSL  routines  (Visual  Numerics,  Inc.,  1984)  were  used  for  random  number  generation. 
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1.3. 1  The  Transition  Number 


The  transition  number,  as  a  sticking  probability,  is  the  only  input  into  SAM.  In 
the  application  of  DLA  to  electrostatics,  Niemeyer,  et  al.  (1984)  and  Wiesmann  and 
Zeller  (1986)  determined  that  the  gradient  of  the  relevant  potential  is  proportionally 
related  to  the  sticking  probability.  In  the  case  of  fluid-fluid  displacement,  we  hypothesize 
that  the  capillary  pressure  is  the  potential  that  we  are  interested  in.  The  gradient  of  the 
capillary  pressure  is  determined  by  combining  the  definition  of  capillary  pressure  with 
Darcy  s  law  written  for  each  phase.  The  resulting  equation  for  the  capillary  pressure 
gradient  is  equivalent  to  the  equation  for  stability  proposed  by  Saffman  and  Taylor  (1958) 
and  Chouke,  et  al.  (1959).  The  interaction  of  the  capillary,  viscous,  and  gravity  forces  at 
the  interface  determines  the  degree  to  which  the  front  is  stable  or  fingered  as  the  NAPL 
displaces  the  water. 

The  capillary  pressure,  p£,  is  defined  as  the  difference  between  the  pressure  in  the 
NAPL  phase  and  the  pressure  in  the  water  phase: 

Pc=PN~Pw  (1-11) 

For  a  two-phase  immiscible  system,  Darcy’s  law  can  be  written  for  each  phase: 


q.  =  -^(,vp.-p.gvz) 

Pw 

(1-12) 

-kkM 

<1n=  (V PN~pNgVz) 

Mn 

(1-13) 

where  w  denotes  water,  N  denotes  NAPL,  and  q  is  the  volumetric  flux  (positive  in  the 
direction  of  gravity).  If  we  differentiate  (1 1)  and  substitute  in  (12)  and  (13),  we  have  the 
equation  for  the  capillary  pressure  gradient. 
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(1-14) 


=}(^-2^)  +  A^sin« 

where  Ap=pN  -  pw ,  and  a  is  the  angle  of  the  flow  from  the  horizontal. 

The  capillary  pressure  gradient  is  then  non-dimensionalized  using  the  same 
definitions  for  non-dimensional  variables  as  Hilfer  and  0ren  (1996).  These  variables  are: 

1  (1-15) 

P  =  PdP  (1-16) 

The  result  of  non-dimensionalizing  the  capillary  pressure  gradient  results  in  Equation  (1- 
17): 


Vp  —  QnPn1  |  AyCg/sina:  (1-17) 

KPa  hPd  pd 

Using  the  previous  definitions  for  the  macroscopic  capillary  and  bond  numbers, 
Equation  (1-17)  can  be  written  as: 


%c=Caw-CaN+Bo  (1-18) 

Morrow  and  Songkran  (1981)  and  Dawson  and  Roberts  (1997)  have  proposed 
similar  linear  combinations  of  the  capillary  and  bond  numbers.  Morrow  and  Songkran 
(1981)  explored  the  effect  of  viscous  and  buoyancy  forces  on  the  entrapment  of  LNAPLs 
in  porous  media.  They  found  that  the  amount  of  entrapped  NAPL  could  be  correlated  by 
a  linear  combination  of  the  microscopic  capillary  and  bond  numbers.  Dawson  and 
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Roberts  (1997),  used  an  approach  similar  to  that  presented  here  to  determine  residual 
DNAPL  saturations,  but  non-dimensionalized  the  capillary  pressure  with  the  Leverett 
(1941)  function: 


V/>c  = 


cos# 

k/0 


VPC 


(1-19) 


In  their  analysis,  the  NAPL  is  assumed  be  immobile. 

The  scaled  capillary  pressure  gradient  in  Equation  (1-18)  is  termed  the  transition 
number.  The  first  two  terms  of  the  transition  number,  T„  compares  the  macroscopic 
capillary  numbers  for  the  wetting  and  non- wetting  phases  and  therefore  captures  the 


combined  effect  of  the  viscous  and  capillary  forces.  The  second  term  in  the  transition 


number,  T2,  is  the  macroscopic  bond  number  and  includes  the  effects  of  gravity  forces 
due  to  the  density  difference  between  the  NAPL  and  water.  The  sum  of  these 
components  is  termed  the  transition  number  as  it  compares  the  relative  importance  of  the 
forces  as  they  transition  between  viscous  dominated,  capillary  dominated,  and  gravity 
dominated.  The  transition  number  is  thus  defined  as: 


T  =  (ca„-CaN)+Bo  (1-20) 

The  macroscopic  capillary  number  has  been  previously  defined.  Bo  is  the  macroscopic 

cw' 

bond  number  (macroscopic  gravity  number  with  the  density  difference  included). 

The  first  term  of  the  equation  is  positive  for  NAPLs  that  are  less  viscous  than 
water  and  negative  for  NAPLs  that  are  more  viscous  than  water.  The  second  term  is 
positive  for  DNAPLs  and  negative  for  LNAPLs,  assuming  the  NAPL  is  displacing  water 
from  the  top.  By  summing  these  components,  the  relative  importance  of  all  three  forces 
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is  taken  into  account.  Positive  values  of  the  transition  number  correspond  to  some  degree 
of  viscous  fingering,  and  negative  numbers  indicate  a  stabilizing  of  the  front.  As  T 
approaches  zero,  the  front  configuration  resembles  a  percolation  cluster. 

1.3.2  The  Transition  Number  as  Model  Input 

The  transition  number  is  the  only  parameter  needed  for  the  model’s  input,  and  is 
based  on  macroscopic  flow  conditions  and  basic  fluid  properties.  It  is  proposed  that  the 
transition  number  determines  the  front  stability  in  the  same  way  that  the  sticking 
probability  changes  the  cluster  density  in  a  SAM  simulation.  This  supposition  is  based 
on  the  work  of  Fernandez  (1991)  and  Kiriakidis  (1991),  but  has  been  defined  differently 
and  has  been  extended  to  include  buoyancy  forces  in  this  work.  As  percolation  clusters 
lie  quantitatively  (in  terms  of  the  fractal  dimension)  and  qualitatively  between  a  fingered 
front  and  a  stable  front,  they  are  correlated  to  a  transition  number  of  0.  The  transition 
number  is  then  mapped  onto  the  range  of  0  to  1,  to  facilitate  its  use  in  the  model  as  a 
sticking  probability,  Ps. 

In  order  to  map  the  ranges  of  T  and  Ps,  it  was  necessary  to  know  the  point  at 
which  the  model  generated  percolation  clusters.  Percolation  clusters  are  characterized  by 
a  fractal  dimension  of  D=91/48  (1.89),  (Berkowitz  and  Balberg,  1993;  Meakin,  1991; 
Oxaal,  et  al.,  1987).  The  box  counting  or  sandbox  method  (Meakin,  1985;  Daccord,  et  al. 
1986),  was  used  to  determine  the  fractal  dimension  of  the  clusters  generated  by  the 
model.  In  this  method,  a  square  of  lattice  sites,  LxL,  is  formed  around  each  occupied  site. 
The  number  of  occupied  sites  within  this  square  is  counted,  N(L).  This  is  done  for 
different  square  sizes  and  a  log-log  plot  of  N(L)  vs.  L  is  generated.  A  straight  line  on  this 
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plot  indicates  that  the  cluster  is  fractal  with  a  fractal  dimension  equal  to  the  slope  of  the 
line.  The  sticking  probability  for  the  model  which  generated  a  percolation  cluster  with 
D=1  89  was  mapped  to  a  zero  value  for  the  transition  number.  This  sticking  probability 
was  determined  to  be  Ps=0.03,  as  the  sticking  probability  changes  as  a  power  law  with 
the  with  cluster  density.  Linear  interpolation  was  then  used  to  relate  the  transition 
number  to  the  sticking  probability  over  its  range. 

1.4  RESULTS  AND  DISCUSSION 

Model  simulations  were  carried  out  for  NAPLs  with  varying  densities  and 
viscosities.  The  angle  of  inclination  of  the  plane  from  the  horizontal  was  90  degrees  and 
15  degrees.  At  each  angle  of  inclination,  two  flow  rates  were  used  to  show  the  effect  of 
the  changing  capillary  and  bond  numbers.  The  NAPLs  used  were  1,2  dichloroethane 
(DC A),  1,2  dibromoethane  (EDB),  gasoline,  and  mineral  oil.  The  first  two  are  DNAPLs 
with  DCA  being  less  viscous  than  water  and  EDB  being  more  viscous  than  water.  The 
last  two  are  LNAPLs  with  gasoline  being  less  viscous  than  water  and  mineral  oil  being 
more  viscous  than  water. 

The  fluid  properties  for  each  NAPL  were  taken  from  Mercer  and  Cohen  (1990) 
and  are  listed  in  Table  1 .  In  Table  2  the  transition  number  as  well  as  its  components,  the 
corresponding  sticking  probability  and  the  fractal  dimension  are  listed  for  each  NAPL  at 
two  different  flow  rates,  and  two  different  angles  of  inclination  from  the  horizontal. 
Figures  5  to  8  show  the  model  simulations  corresponding  to  the  cases  in  Table  2.  In  each 
case,  the  simulation  was  stopped  when  the  finger  or  cluster  reached  the  bottom  of  the 
grid.  Thus  even  though  two  simulations  are  at  the  same  volumetric  flux  rate,  they  may 
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Table  1 

NAPL  properties 


NAPL 

Density  (g/cm3) 

Viscosity  (cP) 

1,2  Dichloroethane 

1.26 

0.887 

1,2  Dibromoethane 

2.17 

1.49 

Gasoline 

0.728 

0.650 

Mineral  Oil 

0.822 

29.0 

Table  2 

Transition  Number  Calculations 


Flux  (cm/s) 

Angle  of  Inclination 

T, 

t2 

T 

Ps 

D 

1,2  Dichloroethane 

0.0076 

90° 

0.0297 

0.3185 

0.3482 

0.3678 

1.67 

15° 

0.0297 

0.0824 

0.1121 

0.1388 

1.81 

0.034 

90° 

0.1331 

0.3185 

0.4516 

0.4681 

1.62 

15° 

0.1331 

0.0824 

0.2155 

0.2391 

1.72 

1,2  Dibromoethane 

0.0076 

90° 

-0.1291 

1.4316 

1.3025 

0.6617 

1.63 

15° 

-0.1291 

0.3705 

0.2414 

0.1471 

1.76 

0.034 

90° 

-0.5774 

1.4316 

0.8542 

0.4443 

1.69 

15° 

-0.5774 

0.3705 

-0.2069 

0.0238 

1.94 

Gasoline 

0.0076 

90° 

0.0922 

-0.3332 

-0.2410 

0.0228 

1.95 

15° 

0.0922 

-0.0862 

0.0060 

0.0358 

1.94 

0.034 

90° 

0.4124 

-0.3332 

0.0792 

0.1068 

1.89 

15° 

0.4124 

-0.0862 

0.3262 

0.3464 

1.68 

Mineral  Oil 

0.0076 

90° 

-7.375 

-0.2181 

-7.593 

0.0235 

1.92 

15° 

-7.375 

-0.0564 

-7.431 

0.0236 

1.94 

0.034 

90° 

-32.991 

-0.2181 

-33.209  0.00153 

2.01* 

15° 

-32.991 

-0.0564 

-33.047  0.00167 

2.01* 

*  D>2  is  a  result  of  the  box  counting  method  used  to  calculate  the  fractal  dimension. 
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not  be  representative  of  the  NAPL  configuration  at  the  same  time.  Work  is  currently 
being  performed  in  which  different  times  for  the  same  flow  rate  are  used  to  show  how  the 
front  progresses. 

These  simulations  show  the  model’s  ability  to  use  one  input  parameter,  the 
transition  number,  to  generate  an  approximate  NAPL  front  configuration.  Instability  of 
the  NAPL  is  shown  in  the  pictures,  as  they  range  from  thin  fingers,  to  fingers  with  thicker 
branches,  to  masses  where  nearly  all  of  the  water  has  been  displaced.  The  configuration 
is  approximate  because,  due  to  random  number  generation,  the  configuration  for  a  given 
transition  number  may  look  somewhat  different  in  each  simulation. 

The  fractal  dimension,  D,  was  used  to  compare  the  relative  density  of  the 
simulations  and  is  included  in  the  tables  for  each  case.  For  1,2  dichloroethane  (Figure  5), 
the  NAPL  configuration  becomes  more  stable,  thicker  branches  in  this  case,  as  the  flow 
rate  is  decreased  and  the  angle  of  inclination  is  decreased.  The  fractal  dimension  ranges 
from  D=1.62  for  the  least  stable  case  at  90°  and  q=0.034  cm/s  to  D=1.81  for  the  most 
stable  case  of  15°  and  q=0.0076  cm/s.  A  flat  front,  where  NAPL  displaces  nearly  all  of 
the  water,  would  have  a  fractal  dimension  approaching  D=2.  The  values  for  the  two 
terms  in  the  transition  number  are  positive  for  each  scenario.  This  indicates  that  higher 
flow  rates  (higher  capillary  number  difference)  and  a  higher  angle  from  the  horizontal 
(higher  bond  number)  both  contribute  to  the  instability  of  a  DNAPL  with  a  viscosity 
lower  than  that  of  water.  The  transition  number  becomes  smaller  as  the  flow  rate  and 
angle  from  the  horizontal  are  reduced,  but  both  terms  remain  positive. 

For  1,2  dibromoethane  (Figure  6),  the  least  stable  case  occurs  when  the  flow  rate 
is  0.0076  cm/s  and  the  angle  of  inclination  is  at  a  maximum.  The  most  stable  case  shown 
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here  results  when  the  flow  is  increased,  but  the  angle  is  decreased.  Here,  the  fractal 
dimensions  have  a  larger  range  spanning  from  D=1.63  to  D=1.94.  When  evaluating  the 
transition  number  components,  we  see  this  time  that  the  viscous  and  gravity  forces  are  in 
opposition  due  to  the  relative  viscosity  of  EDB  to  water.  For  larger  flow  rates,  the  first 
term  in  the  transition  number  becomes  more  negative  indicating  that  for  a  DNAPL  with  a 
viscosity  greater  than  that  of  water,  higher  flow  rates  will  contribute  to  the  fronts 
stability.  For  larger  angles  of  inclination,  the  second  term  in  the  transition  number 
becomes  greater  indicating  more  unstable  flow.  For  this  NAPL,  the  forces  oppose  each 
other. 

The  third  NAPL,  gasoline  (Figure  7),  shows  a  trend  opposite  that  of  EDB. 

Because  gasoline  is  less  viscous  than  water,  it  becomes  increasingly  unstable  as  the  flow 
rate  is  increased.  Stability  is  increased  in  this  case  by  increasing  the  angle  of  inclination 
due  to  the  negative  density  difference  between  the  two  fluids.  A  crossover  between 
inherently  stable  and  inherently  unstable  front  configurations  can  occur  when  the  terms  of 
the  transition  number  add  to  zero.  This  is  approached  when  q=0.0076  cm/s  and  the  angle 
from  the  horizontal  is  15°. 

In  the  final  case,  mineral  oil  (Figure  8)  is  used  to  show  the  front  configuration  for 
a  NAPL  which  is  less  dense  but  more  viscous  than  water.  For  this  scenario,  both  terms  in 
the  transition  number  are  negative  regardless  of  the  changing  flow  rate  and  angle  from  the 
horizontal.  The  fractal  dimensions  of  the  simulations  range  from  D=1.92  to  D=2.01. 

With  a  viscosity  29  times  greater  than  that  of  water,  the  first  term  of  the  transition  number 
dominates  the  displacement  behavior  and  the  differences  due  to  the  angle  of  inclination 
are  small.  A  fractal  dimension  greater  than  two  is  a  consequence  of  the  inaccuracy  of  the 
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box  counting  method  used  to  calculate  the  fractal  dimension.  Specifically,  when  the 
cluster  reaches  the  sides  of  the  grid,  the  box  counting  method  does  not  account  for  the 
filled  grid  spaces  at  the  very  edge  of  the  cluster. 

The  simulations  for  each  NAPL  show  approximate  configurations  for  NAPL- 
water  displacement  in  a  porous  medium.  This  configuration  is  a  result  of  evaluating  the 
viscous,  gravity  and  capillary  forces  at  the  pore  level.  The  simulations  presented  here  do 
not  show  a  residual  saturation,  but  rather  a  percentage  of  pore  space  in  which  the  NAPL 
displaced  the  water.  Using  this  information,  retention  capacity  values  and  subsequent 
penetration  depths  can  be  determined  more  accurately.  Future  work  for  this  model 

includes  the  incorporation  of  simple  heterogenieties,  and  the  calculation  of  retention 
capacities. 
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Figure  5.  The  displacement  of  water  by  1,2  Dichloroethane  at  different  flow  rates  and  different 
angles  of  inclination  from  the  horizontal:  (a)  q=0.0076  cm/s,  <x=90°,  (b)  q=0.034  cm/s,  a=90°, 
(c)  q=0.0076  cm/s,  a=15°,  (d)  q=0.034  cm/s,  a=15°. 
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(c)  q-0.0076  cm/s,  a=l 5°,  (d)  q=0.034  cm/s,  a=l 5°. 
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Figure  7.  The  displacement  of  water  by  gasoline  at  different  flow  rates  and  different  angles  of 
inclination  from  the  horizontal:  (a)  q=0.0076  cm/s,  a=90°,  (b)  q=0.034  cm/s,  a=90°,  (c) 
q=0.0076  cm/s,  a=15°,  (d)  q=0.034  cm/s,  a=T5°. 
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Figure  8  .  The  displacement  of  water  by  mineral  oil  at  different  flow  rates  and  different  angles 
of  inclination  from  the  horizontal:  (a)  q=0.0076  cm/s,  <x=90°,  (b)  q=0.034  cm/s,  <x=90°,  (c) 
q=0.0076  cm/s,  a=15°,  (d)  q=0.034  cm/s,  a=l  5°. 
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1.5  CONCLUSIONS 

A  modeling  technique  for  describing  the  front  configuration  for  NAPLs 
displacing  water  in  saturated  porous  media  is  presented.  The  stochastic  aggregation 
model  (SAM)  is  a  modified  diffusion  limited  aggregation  model  that  can  quickly  provide 
two-dimensional  representations. 

Capillary,  viscous  and  gravity  forces  act  concurrently  to  stabilize  or  destabilize 
the  front  when  a  NAPL  displaces  water.  A  transition  number  has  been  derived  that 
determines  the  interplay  of  these  forces,  and  is  a  scaled  linear  combination  of  the 
macroscopic  capillary  and  bond  numbers.  The  use  of  a  transition  number  as  the  model 
input  is  unique  in  that  it  includes  gravity  as  well  as  the  fluid’s  viscosity  and  density 
differences  and  allows  the  model  to  create  simulations  ranging  from  highly  fingered  to 
stable  displacements.  Easily  measured  or  estimated  macroscopic  parameters  are  used  to 
calculate  the  transition  number,  thereby  contributing  to  the  ease  of  use  and  broad 
application  of  the  model. 

In  this  study,  we  presented  model  simulations  for  four  different  NAPLs:  two 
LNAPLs,  with  one  less  and  one  more  viscous  than  water,  and  two  DNAPLs,  with  one 
less  and  one  more  viscous  than  water.  The  effects  of  different  flow  rates  and  different 
inclinations  from  the  horizontal  were  investigated  for  each  NAPL.  Fractal  dimensions 
were  calculated  and  used  to  compare  model  results.  The  simulations  accurately  show  that: 

•  for  a  DNAPL  which  is  less  viscous  than  water,  as  both  the  NAPL  flow  rate  and 
the  angle  from  the  horizontal  is  increased  the  front  becomes  less  stable; 
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•for  a  DNAPL  which  is  more  viscous  than  water,  as  the  NAPL  flow  rate  increases 
the  front  becomes  more  stable,  but  as  the  angle  from  the  horizontal  increases  the  front 
becomes  less  stable; 

•for  a  LNAPL  which  is  less  viscous  than  water,  as  the  NAPL  flow  rate  increases 
the  front  becomes  less  stable,  but  as  the  angle  from  the  horizontal  increases  the  front 
becomes  more  stable; 

•for  a  LNAPL  which  is  more  viscous  than  water,  as  both  the  NAPL  flow  rate  and 
the  angle  of  inclination  is  increased  the  front  becomes  more  stable. 
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SECTION  V 


SECTION  V 


SUMMARY  AND  CONCLUSIONS 

1.1  INTRODUCTION 

This  project  focused  on  the  distribution  of  multiphase  fluids  in  homogeneous  and 
layered  porous  media  under  equilibrium  and  transient  flow  processes.  A  synergistic 

procedure  was  used  that  included  pore-scale  and  lab-scale  numerical  models,  and  one-  and 
two-dimensional  laboratory  studies. 

1.2  CONCLUSIONS 

This  report  is  divided  into  several  sections.  Section  I  is  the  introduction.  Section  II 
reports  results  from  the  physical  one-and  two  dimensional  laboratory  studies.  Section  m 
presents  and  verifies  a  one-dimensional  transient  model  for  multiphase  flow  through 
layered  porous  media.  Finally,  Section  IV  presents  a  pore-scale  modified  diffusion  limited 

aggregation  model.  Conclusions  are  listed  at  the  end  of  each  chapter.  Major  conclusions 
are  summarized  here. 

1)  Effects  of  soil  layering,  nonwetting  phase  access  and  differences  in  those  effects 
for  air  and  LNAPL  systems  were  quantified  in  a  one-dimensional,  transient 
drainage  experiment  in  layered  soil. 

2)  Perfluorocarbons  have  fluid  properties  that  are  similar  to  those  of  dense, 
chlorinated  solvents  but  pose  no  health  hazard  and  are  not  regulated.  These  non- 
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toxic,  inert,  dense  nonaqueous  phase  liquids  were  identified  and  investigated  as 
surrogates  for  laboratory  and  field  studies  of  chlorinated  solvents. 

3)  Stable  mounds  of  DNAPL  can  form  on  heterogeneities,  with  lateral  stability 
existing  on  flat,  or  even  sloped,  textural  discontinuities. 

•  A  predictive  expression  was  developed  for  the  maximum  thickness 
of  a  laterally  stable  NAPL  mound  and  it  was  found  that  hysteretic 
capillary  pressure-saturation  relations  are  needed  to  explain  LNAPL 
and  DNAPL  mound  thicknesses  greater  than  zero. 

•  Models  which  do  not  include  hysteresis  will  not  result  in  correct 
distributions  of  free  phase  NAPL  distribution  at  equilibrium. 

4)  A  conceptual  model  of  LNAPL  infiltration  and  spreading  at  the  capillary  fringe 
was  developed  and  verified  with  a  series  of  two-dimensional  flume  experiments. 

•  The  top  of  the  lens  is  not  necessarily  at  the  top  of  the  capillary  fringe  but  its 
location  can  be  calculated  from  the  water-NAPL  and  NAPL-air  interfacial 
tensions. 

•  Capillary  pressures  at  each  point  on  the  water-NAPL  boundary  are  on  a 
different  hysteresis  loop  of  the  water-NAPL  capillary  pressure-saturation 
curve. 

•  The  top  of  an  LNAPL  lens  is  flat  and  the  capillary  pressure  at  each  point  is 

the  air-entry  pressure  on  the  NAPL-air  drainage  curve. 

•  The  maximum  thickness  of  an  LNAPL  lens  is  the  difference  between  the 
capillary  pressure  entry  head  on  the  oil-water  drainage  curve  and  the 
capillary  pressure  entry  head  on  the  oil-water  imbibition  curve. 

5)  Both  a  two-dimensional  modified  diffusion  limited  model  and  a  three-dimensional 

invasion  percolation  model  were  developed  as  part  of  this  project.  Results  include: 

The  modified  diffusion  limited  aggregation  model  appears  to  simulate  the 
unstable  flow  of  DNAPLs  in  saturated  porous  media  well.  We  found  that 
the  model  was  able  to  simulate  both  DNAPL  and  LNAPL  in  two- 
dimensions. 
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•  The  volume  of  the  entrapped  nonwetting  phase  is  determined  mainly  by  the 
aspect  ratio  of  the  network  and  pore  tortuosity.  As  the  aspect  ratio 
decreases,  entrapment  of  the  nonwetting  phase  and  hysteresis  decreases. 

•  The  bond  number  was  found  to  play  a  major  role  in  determining  the 
drainage  and  imbibition  sequences  and  hence  strongly  influence  the  water 
retention  and  hydraulic  conductivity  curves.  This  study  showed  that,  as  the 
bond  number  increased,  the  relative  hydraulic  conductivity  and  water 
retention  curves  become  scale  dependent  for  coarse  soils. 

6)  A  semi-analytical  solution  for  one-dimensional,  two-fluid  flow  in  layered  porous 
media  is  presented.  The  modeling  approach  presented  in  this  report  has  the 
following  properties:  a  set  of  equations  is  solved  explicitly  and  sequentially,  rather 
than  simultaneously  as  required  by  the  finite  difference  and  finite  element  methods, 
spacing  between  nodes  is  varied  within  and  between  time  steps  based  on  accuracy 
criteria,  and  the  iteration  process  requires  convergence  of  only  one  or  two 
boundary  condition  values,  rather  than  iterations  on  all  internal  nodes. 

•  The  model  was  found  to  be  unconditionally  stable  for  one-dimensional, 
two-fluid  transient  flow  in  layered  soil. 

•  The  model  was  verified  by  comparing  model  results  to  an  exact  solution 
and  by  comparing  simulations  to  laboratory  data. 

•  A  comparison  of  Richards’  equation  and  two-fluid  flow  equations  for 
drainage  with  limited  air  access  showed  that  the  hydraulics  of  air  flow  can 
have  a  significant  impact  on  the  nature  of  vertical  water  drainage. 
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